Remote identification of non-lambertian materials

ABSTRACT

In one example of a method for remote identifying a non-Lambertian target material, a spectral signature for a target is determined from each of at least two different sets of imagery acquired at different angles, and compared to a predicted signature for a candidate material for each of the at least two different angles. The predicted signatures take into account the known anisotropy of reflectance, and thus also radiance, of the candidate material.

CLAIM OF PRIORITY

This patent application claims the benefit of priority of U.S. patentapplication Ser. No. 12/851,991, filed on Aug. 6, 2010, which is herebyincorporated by reference herein in its entirety.

TECHNICAL FIELD OF THE INVENTION

The invention relates generally to remote sensing and materialidentification.

BACKGROUND OF THE INVENTION

Material identification using a remotely located sensor is used for anumber of purposes, for example detection, identification andclassification of objects within a scene, and applications, forcharacterization of urban areas, search and rescue, differentiatingcombatant forces, and detection of attempts at camouflage, denial anddeception. It is based on the general principle that the observedspectral radiance of any given target will vary based on, among otherthings, the type of the material of which the surface of the target ismade. Different materials absorb, reflect and emit differently,depending on wavelength. The spectral radiance of a target, the targetbeing the surface of a particular object within a scene, from a givenangle or direction can be measured using various types of sensors,depending on the wavelengths of interest. Based on the measured orobserved spectral radiance it is possible to determine the material ofwhich the surface is made.

Several types of remote sensors and imaging modalities have been used togenerate image sets containing both spectral and spatial information ofreal scenes for purposes of detection, identification or classificationof objects within the scene. Electro-optical sensors are typicallysingle band, multispectral or hyperspectral. A multispectral sensordetects and records radiation in a limited number of bands that aretypically fairly wide, for example in the red, blue, green, andnear-infrared bands of the spectrum. A hyperspectral sensor detectsradiation in a large number of contiguous bands, typically throughoutthe visible and near-infrared regions of the spectrum. Other types ofimaging modalities, for example, synthetic aperture radar (SAR),typically operate only in a single band. The sensors are typically (butdo not have to be) placed in satellites or aircraft and acquire imagesof portions of the surface of the earth during flyovers at relativelyhigh altitudes. However, it is possible for the sensors to be placed onthe ground.

Each “image”—also called “imagery” or “image set”—of a scene generatedby such a sensor comprises spatial information, typically in twodimensions. It also contains spectral radiance information, which wouldinclude the radiance of at least one predetermined band of wavelengthsthat the sensor can detect. The material of which at least the surfaceof an object within the scene is made, called the “target,” isidentified by selecting within the image pixels comprising the targetand then evaluating the spectral radiance of those pixels to develop aspectral signature for the target that can be compared to spectralsignatures for various materials. In automatic material identification,a specially programmed computer is used to process image data from aremote sensor and other data to identify the material of the target.

The spectral radiance, which is radiance at a given wavelength or band,for any given target in a scene will depend on the material of which thetarget is made, as well as the spectrum and angle of irradiation beingreflected by the target, the atmospheric conditions through which boththe illuminating irradiation and the reflected radiation travels, andthe spectrum of any emitted radiation. In order to make theidentification, measured spectral radiance is typically transformed toan estimated reflectance. Reflectance is the ratio of the measuredradiance from an object divided by the radiance reflected by a 100%Lambertian reflector. When using images of real targets, reflectance issometimes estimated by taking into account relevant environmentalconditions, such as the radiation source and atmosphere, under which theimagery was acquired.

The way in which a surface reflects or emits radiation can be generallycategorized as either Lambertian or non-Lambertian. A Lambertian surfacescatters electromagnetic radiation equally in all directions, withoutregard to the direction of illumination. Thus, its reflectance isgenerally isotropic, or the same in all directions. A non-Lambertiansurface does not scatter incident electromagnetic radiation equally inall directions. Examples of non-Lambertian surfaces include those thatare backscattering, meaning that the light scatters predominantly towardthe illumination source; forward scattering, meaning scatteringpredominantly in directions away from the illumination source; andspecular, meaning reflecting the illumination source like a mirror. Manyman-made objects or targets exhibit non-Lambertian reflectance.

In order to identify a material within an image, prior art methods treatthe material as Lambertian. The candidate materials are also treated asLambertian. The directional hemispherical reflectance (DHR) for eachcandidate material, which relates, for a given wavelength or band ofwavelengths and direction of incident irradiation, reflected radianceacross the entire hemisphere, is used to predict the spectralreflectance of the candidate material.

SUMMARY

Treating all targets as Lambertian in remote material identificationprocesses tends to yield erroneous results or no results at all whenremotely identifying non-Lambertian materials using spectral signatures.The invention generally relates to methods and apparatus useful inconnection with identification using remote sensors of non-Lambertiantarget materials.

According to one example of a method for remotely identifying anon-Lambertian target material, a spectral signature for a target isdetermined from each of at least two different sets of imagery acquiredat different angles, and compared to a predicted signature for acandidate material for each of the at least two different angles. Thepredicted signatures take into account the known anisotropy ofreflectance and/or emissivity of the candidate material. If the targetmaterial is non-Lambertian, using two angles and predicted spectralsignatures for candidates for those two angles tends to increase theprobability of correct identification of the target material and tendsto decrease the probability of false matches.

According to another example of a method for remotely identifyingmaterials, predicted spectral signatures for non-polarimetric,reflective imagery are determined from bi-directional reflectancedistribution functions and do not rely solely on DHR.

According to yet a different example of a method for remote identifyingmaterials, an analysis of uncertainties in an estimated multi-angletarget spectral signature and a multi-angle candidate spectral signatureis used to determine level of confidence in a decision on a match.Correlated errors in the multi-angle spectral signatures of the targetand the candidate material are, preferably, considered.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1A is a schematic illustration of an area containing a target ofinterest, being illuminated by the sun and imaged by two remote, passivesensors.

FIG. 1B is a schematic illustration of the area of FIG. 1A being imagedby two remote, active sensors.

FIG. 1C is a schematic illustration of the area of FIG. 1A with anemissive target being imaged by two sensors.

FIG. 2 is a flow diagram of a method for identifying a material usingremote sensors.

FIG. 3 is a schematic illustration of a representative computing systemspecially programmed to execute the methods for identifying materialsusing remote sensors according to the method of FIG. 1, and any one ormore of examples of methods for different imaging modalities describedbelow.

FIG. 4 is a schematic illustration of one representative implementationof programs executing on the computing system of FIG. 3 when it isexecuting certain steps of the process of FIG. 2.

FIGS. 5A and 5B are graphs illustrating a graphical comparison ofmulti-angle target signatures with multi-angle candidate signatures.

FIG. 6 is an example of an embodiment of a remote materialidentification process for non-polarimetric reflective imagery.

FIG. 7 is a second example of an embodiment of a remote materialidentification process for non-polarimetric reflective imagery.

FIG. 8 is an example of an embodiment of a remote materialidentification process for non-polarimetric emissive imagery.

FIG. 9 is an example of a process for calculating a best-fittingcandidate spectral radiance signature for the process of FIG. 8.

FIG. 10 is second example of a process for calculating a best-fittingcandidate spectral radiance signature for the process of FIG. 8.

FIG. 11 is an example of an embodiment of a remote materialidentification process for polarimetric reflective imagery

FIG. 12 is an example of an embodiment of a remote materialidentification process for polarimetric emissive imagery.

FIG. 13 is an example of an embodiment of a remote materialidentification process for synthetic aperture array imagery.

DETAILED DESCRIPTION

In the following description, like numbers refer to like elements.

Referring to FIGS. 1A, 1B and 1C, each figure schematically illustratesdifferent examples of acquiring multi-angle imagery. Sensors 102 and 104each acquire imagery of at least a portion of an area of earth 106—ascene—containing one or more targets of interest, for example, target108. In the example of FIG. 1A, the sun 110 is illuminating the scene.However, other sources of illumination of the target are possible,including, for example, other man-made sources of illumination. FIG. 1Billustrates the case of the sensors 102 and 104 being active sensors,meaning that they illuminate the scene with radiation and measurereflectance of that radiation. Examples include synthetic aperture radar(SAR) systems and light detection and ranging (LIDAR) systems. An activesensor, such as used for mono-static SAR, is illustrated by FIG. 1B.However, the illumination source need not necessarily be co-located withthe sensor, as in multi-static SAR. Radiance from a target may alsoinclude emitted radiation, in particular radiation in the infraredwavelengths, which is dependent on the target's temperature. Acquisitionof imagery by sensors for detecting such emitted radiation isillustrated by FIG. 1C.

Aircraft or satellites, for example, carry the sensors. However, thesensors could be placed anywhere that affords a view of the areacontaining a target of interest. Images acquired from the sensors arestored on computer readable media for further processing. The data istransmitted or delivered to an image processing facility (not shown) forreconstruction (if necessary) and further processing. The same or adifferent facility can perform material identification using theprocessed imagery.

Sensors 102 and 104 that are positioned in each of FIGS. 1A, 1B and 1Care intended to represent two different image acquisition angles. Theacquisition angles, also referred to as the acquisition geometry are,unless the context otherwise indicates, angles in three dimensions thatare typically defined by a zenith angle (measured from the targetsurface normal) and an azimuthal angle. Different image acquisitionangles can differ in terms of zenith angle, azimuthal angle, or both.Although they are illustrated as two separate sensors, sensors 102 and104 in each of the figures are also intended to represent the situationin which the same sensor moves positions. They can also be differentsensors of the same type or different types of sensors.

In the example that is illustrated in FIG. 1A, sensors 102 and 104 arepassive electro-optical sensors. The sensors can be single band,multispectral, or hyperspectral. Sensor 102 acquires imagery of thetarget from along a first direction 116, and the second sensor 104acquires an image of the target from a second direction 120. “Direction”refers to the direction of propagation of radiation from the target tothe sensor or, in the case of the illumination source, from thedirection of propagation of irradiation from the source to the target.In the case of sensor 102, the zenith angle 114 is measured between thenormal 112 and direction 116, and, in the case of sensor 104, it is theangle 118 measured between normal 112 and direction 120. The azimuthalangle for sensor 102 is angle 122, and for sensor 104, it is angle 124.In this example, the plane of the surface of the target is presumed tobe the plane of the terrain 106, but it does not need to be.

Referring now also to FIG. 2, flow chart 200 illustrates the basic stepsof a method for identifying the material, such as target 108 of FIG. 1,using images taken from two or more acquisition geometries. Certain ofthe steps are performed by a specially programmed computing system. Theprocessing steps are the same for both single-modality imagery (imagesusing the same modality) and multi-modality (images taken usingdifferent modalities) operation.

As represented by step 202, two or more images containing the sametarget are acquired at step 202 using remote imaging sensors. The imagesneed to each contain the target, but otherwise do not need to becoextensive, or to be taken using the same type of sensor or imagingmodality. At least two of the images are acquired from different angles.The sensors can be, for example, mounted in aircraft or satellitesmoving across the area or, in some situations, positioned on the ground.

Also obtained is the information about the environmental conditionsunder which the images were acquired. The particular information that isobtained depends on the imaging modality. Environmental conditionstypically include atmospheric conditions, which can include aerosol,water vapor content, and temperature. Environmental information can beacquired in connection with the acquisition of each image. It can alsobe obtained subsequent to the acquisition of the images usingcontemporaneous observations recorded in other databases. The particularenvironmental information that is acquired depends at least in part onthe modality of the imagery.

More than two sensors can be used to acquire images of a target. Thesame sensor can also be used to acquire images from each of the angles.If different sensors are utilized to acquire the two or more images ofthe target, the sensors need not utilize the same imaging modality. Theacquired images are typically recorded on local media. The images arecommunicated to a remote facility for processing and use by transportingthem on computer readable physical media, transmitting them usingsignals, or by combination of communication modes.

In the example of FIG. 1A, sensors 102 and 104 are of a type that senseradiation in the visible and near infrared bands that has been reflectedby the target material. The source of irradiation incident on the targetis, in this example, the sun 110, but there could be artificial sourcesilluminating the target. In this example, sensor 102 is positioned in a“near specular” geometry and sensor 104 is positioned in an “offspecular” geometry. Sensor 102 is positioned to detect the spectralradiance within the forward scattering lobe of the radiation from sun110 reflected by target 108. The width of the forward scattering lobe,and thus the range of angles (azimuthal and zenith) at which the sensoris positioned, depends on the particular material of which the surfaceof the target is made. The sensor 102 need not be located exactly in thespecular direction, which is 180 degrees in azimuth from the sun and thesame zenith angle as the direction of irradiation on the target. Infact, depending on the type of sensor and illuminating source, placing asensor in this position could result in the sensor being saturated andthus unable to provide spectral information. Sensor 104 is acquiring animage that is used to determine spectral radiance at an acquisitionangle outside the forward scattering lobe. Either of the two sensorscould be, instead, located within the back scattering lobe of reflectedradiation. Spectral radiance measurements within the forward scatteringlobe and outside the forward scattering lobe will tend to give the bestidentification results for reflected radiation of non-Lambertiantargets.

For active imaging modalities in which the illuminating source and thesensor are co-located, such as mono-static SAR and LIDAR, betteridentification results for non-Lambertian materials tend to be achievedwhen at least one of at least two images are, for example, taken nearnormal orientation to the target surface and at least one is taken at anelevation off the normal orientation. This is illustrated by FIG. 1B.Multi-static SAR, in which the source of the microwave radiation is notco-located with the sensor, and images acquired of artificiallyilluminated scenes are treated like the reflective imaging modality,such as shown in FIG. 1A. The microwave radiation source would betreated similarly to the sun 110. For emissive imaging modalities, suchas infrared, in which the target is emitting radiation, taking one imagefrom a direction that is near normal to the surface, and another takenat an elevation far from normal to the surface tends to give betteridentification results for non-Lambertian targets. This is illustratedby FIG. 1C.

The target is delineated in each of two or more images at step 204 byselecting the pixels in each image that comprise a target of interest.The selection can be done manually or automatically by a computingsystem based on predetermined criteria.

At step 206 the spectral radiance of the target, which can be comprisedof both reflected and emitted radiation, is obtained by a speciallyprogrammed computing system from at least two of the images that areacquired at different angles. The spectral signature of the target anderror bars for the signature are estimated for each image at step 208 bythe specially programmed computing system. The errors bars are based onthe errors sources associated with the target signature. The targetsignature for each image is then concatenated into a multi-angle targetspectral signature.

Table 1 lists the different imaging modalities and, for each modality,the type of target signature that is estimated, and sources of error inthe determining the estimated target signature and its error bars.

TABLE 1 Imaging Modality Target Signature Error Sources Non-polarimetricReflectance or Atmospheric conditions (non-PI) reflective radiancespectrum Radiometric calibration single band or BRDF (bi-directionalmultispectral (MS) reflectance distribution function) measurementsNon-PI reflective Reflectance or Same as non-PI reflective hyperspectral(HS) radiance spectrum MS, plus: Spectral calibration Non-PI emissiveRadiance spectrum Same as non-PI reflective MS single-band or MS Non-PIemissive Radiance spectrum Same as non-PI reflective HS MS, plus:Spectral calibration Polarimetric (PI) Radiance spectrum Same as non-PIreflective reflective or MS, plus: emissive single- PI calibration bandor MS PI reflective or Radiance spectrum Same as PI reflective emissiveHS MS, plus: Spectral calibration Synthetic Aperture Reflectance or Sameas non-PI reflective MS Radar (SAR) radiance spectrum

Although the type of target signature for non-polarimetric reflectiveimagery is indicated to be spectral reflectance, the type of signatureused in the process can either be spectral radiance or spectralreflectance. Spectral reflectance signatures can be determined from theimage radiance in two ways. If the sensor output is linearly related toradiance and the scene contains horizontal Lambertian objects of knownreflectance, the empirical line method can be used. The other way is touse model-based atmospheric compensation terms derived from atmosphericparameters and background reflectance.

One or more candidate materials are selected at step 210. For eachcandidate material, as indicated by step 212, a candidate spectralsignature is estimated at step 214 by the specially programmed computingsystem for each image based on acquisition geometry, the imagingmodality, surface reflection parameters from a bi-directionalreflectance distribution function (BRDF) database for candidatematerials, and relevant imaging conditions, which include background andatmospheric information. The candidate spectral signatures areconcatenated into a multi-angle candidate spectral signature. If theestimated target spectral signature for particular acquisition geometryis determined in the radiance domain, the corresponding predictedcandidate spectral signature is also determined and expressed in theradiance domain. Similarly, for estimated target spectral signaturesexpressed in reflectance, the predicted candidate spectral signaturesare determined in the reflectance domain.

Candidate spectral signatures are determined by taking into account thereflectance anisotropy, which is defined for a given material by theBRDF. If the imaging modality is reflective, the reflectance signatureof the candidate material is calculated from its Lambertian-equivalentsolar reflectance signature and directional hemispherical reflectancesignature from the BRDF database. BRDF is the ratio of reflectedradiance to incident irradiance at a given wavelength, and at givendirections of propagation. The direction is characterized by a zenithangle with respect to the normal to the surface of the material, and anazimuth angle, as indicated in FIG. 1. A BRDF for a given material canbe empirically measured for all wavelengths and directions, or it can bepredicted from a BRDF model for the material, and thus has errorassociated with it.

The predicted candidate signatures are then used in a decision processto decide if the candidate material matches the target within apredetermined level of confidence. At step 216 the specially programmedcomputing system calculates the differences between the estimated targetmulti-angle signature and the predicted candidate spectral signature.Steps 214 and 216 are, as indicated by steps 212 and 218, repeated foreach candidate material. A decision is made at step 220 as to which, ifany, of the candidate materials match the target material to within thepredetermined level of confidence. Many approaches to the decisionprocess are possible. One example is a chi-square test. When a candidatesignature matches the target signature closely enough, the candidatematerial is determined by the process to be a possible target materialand is listed at step 220 as such by the computing system executing theprocess.

FIG. 3 illustrates a representative, specially programmed computingsystem 300 for executing software that performs at least the functionsdescribed herein in connection with steps 206, 208, 210, 214, 216, 218,and 220 of the material identification process of FIG. 2. The computingsystem 300 comprises one or more processors, which are represented byprocessing entity 302, and working memory 304, in which programinstructions are loaded for execution by processing entity 302. Theseinstructions can also be stored or contained on any type of electronic,optical or magnetic media, as well as on transient signals, that can beread by the computing system. The instructions, which may first need tobe compiled, are loaded into the memory 304 for execution by theprocessor. Examples of such media include mass data storage andremovable media. Information and data on which the process acts, as wellas resulting from operations of the processing entity 302, are stored inmemory 304.

The system also includes an input/output subsystem 306, which isrepresentative of one or more subsystems through which the computingsystem may interact with a user or may communicate with other computingsystems by transmitting information using signals. Examples of the oneor more subsystems include a display, a user input device, such as akeyboard, mouse, touch pad, touch screen, or remote gesture recognitiondevice, through which a user may interact with the program, andinterfaces for communicating with other computing systems or devices. Noparticular computer architecture is intended to be implied by thisexample. The example is intended to be representative generally ofcomputing systems suitable for being programmed to perform theseprocesses, and not limiting. Furthermore, the processing need not belimited to a single computing system, but could be distributed amongmore than one computing system. Furthermore, different programs runningon computing system 300 or on multiple computing systems may executeparts of the process described by FIG. 2.

The process carried out the computing system 300 operates on imagerystored in imagery storage system 308. Storage system 308 is intended tobe representative of any type of system for storing acquired imagery. Itcould be a simple file system or one or more databases, for example.Imagery acquired through, for example, sensor 102 and 104 of FIG. 1would be stored in system 308. The storage system 308 is illustrated asbeing local, communicating over bus 310 with the processor 302 andmemory 304. It could be stored on a local hard drive, for example.However, it could also be located on a mass data storage device locatedelsewhere on a local or wide area network, or on a remote server. BRDFdatabase 312 is representative of one or more databases containingpredetermined or previously measured information or models fordetermining direction and wavelength dependent surface reflectancecharacteristics and parameters of candidate materials. It is, likeacquired imagery 308, shown as being stored locally; it communicatesover bus 310 with the processor. However, it can be stored on a local orwide area network, or accessed remotely through a server.

Referring to FIG. 4, illustrated is an example of one possiblearrangement of programs that can be used to implement the process ofFIG. 2 and the processes described below on a computing system, such ascomputing system 300. Program 402 carries out the materialidentification process generally illustrated by FIG. 2. However, itmakes use of, among others, two modeling programs to carry out severalof the steps. Program 404 is called MODTRAN®, which is a program ownedby the U.S. Air Force Research Laboratory, maintained by SpectralSciences Inc. and distributed by Ontar Corporation. It modelsatmospheric propagation of electromagnetic radiation, including theeffects of molecular and particulate absorption/emission and scattering,surface reflections and emission, solar/lunar illumination, andspherical refraction. MODTRAN® is used by the basic program to determinehow the atmosphere affects the radiation incident on, and reflectedfrom, the target. Program 406 is a suite of publicly available programscalled the Nonconventional Exploitation Factor Data System (NEFDS). TheNonconventional Exploitation Factor (NEF) database 408 is a database ofBRDF for a large number of materials. Information about the NEFDSprograms and the NEF database was developed by the Cortana Corporationfor the National Geospatial-Intelligence Agency (NGA) of the UnitedStates. The NEFDS programs use the database to calculate, among otherthings, BRDF, while taking into consideration the properties of thetarget material, the material background, the measuring sensors, and theatmospheric environment. The database contains pre-measured surfacereflection parameters for over 400 materials corresponding to a widevariety of objects ranging from camouflage to paint. A NEFDS functioncan be used to query for BRDF values of any of these materials at arange of wavelengths and any given geometry, and is used by the processfor, among other things, predicting candidate spectral signatures andtransforming target radiance signatures into target reflectancesignatures. No particular programming arrangement, architecture orframework is intended to be implied by the foregoing examples.

FIGS. 5A and 5B are graphs with reflectance on the vertical axis andspectral bands on the horizontal axis displaying representations of aconcatenated, multi-angle spectral signature of a target and of severalcandidate materials. FIG. 5A illustrates two examples where incorrectcandidate material signatures are compared to the target signature whichis represented in each figure by error bars 502 (off-specular image) and504 (near-specular image). The signature of the first incorrectcandidate is represented by error bars 506 for the off-specular imageand error bars 508 for the near specular image. If only thenear-specular image was used for identification matching, the error bars506 are close enough that the first incorrect candidate material wouldresult in a false match. However, when both specular and off-specularimages are considered, the differences between the signatures in theoff-specular image prevent the target from being confused with the firstincorrect candidate material. The signature of the second incorrectcandidate is represented by the error bars 510 (off-specular) and 512(near specular). In this case, the signature of the second incorrectcandidate matches the target signature in the off-specular image but notin the near-specular image. The use of both images again suppresses thesecond incorrect material as a false match. FIG. 5B graphicallyillustrate a match between a candidate material's signature, representedby error bars 514 (off-specular) and 516 (near-specular).

Following are descriptions of examples of implementations of the processof FIG. 2 for different imaging modalities.

Non-Polarimetric Reflective Imagery

FIGS. 6 and 7 illustrate two examples of implementation of the processof FIG. 2 for non-polarimetric reflective imagery which are performed onone or more specially programmed computers. The reflective imagery maybe a single band, multispectral or hyperspectral. These examples assumemultispectral imagery. Non-polarimetric hyperspectral imagery isprocessed in a substantially similar manner. FIG. 6 is an example ofperforming the process in the reflectance domain, and FIG. 7 is anexample of performing the process in the radiance domain.

Each process uses target radiance from at least two images. The imagesare generated using at least one radiometrically calibratednon-polarimetric multispectral sensor. The images are acquired at two ormore acquisition geometries, as generally shown in FIG. 1A. In thefollowing example, at least one image is acquired near the speculargeometry and one acquired far from the specular geometry, and the targetis assumed to be a non-Lambertian material.

For each image, each of these examples utilizes the following inputs:target radiance; background radiance; the relative spectral responsefunction for each band in the images data set; information about theatmosphere through which the target was illuminated by radiation andthrough which the reflected radiation was sensed; candidate materialsurface properties; and uncertainties.

Information about the atmosphere includes the effective aerosolparameters for each aerosol regime, primarily the boundary layer, thetroposphere and the stratosphere, as well as other atmosphericparameters and their uncertainty, for example, water vapor columnabundance. Examples of the aerosol parameters include effectiveextinction, absorption and asymmetry parameters for an assumed layerthickness, and an associated error covariance matrix. One or more of theparameters can be used. Total column amounts are usually sufficient;vertical profiles are not necessary. The background radiance is used toderive background reflectance. Candidate material surface propertiesinclude diffuse reflectance and Lambertian-equivalent solar reflectance.They can be obtained from the BRDF database 602, such as through the NEFdata system. Uncertainty inputs include radiance measurementuncertainty, error covariance matrix for the diffuse reflectance andLambertian-equivalent solar reflectance, aerosol parameter uncertainty,and spectral calibration uncertainty for hyperspectral sensor imagery.

The following steps are performed by one or more specially programmedcomputing systems for each image that is being used in theidentification process. Aerosol parameters 604 are used at step 606 toestimate background independent atmospheric correction for each image.The aerosol parameters are also used to estimate a “BRDF atmosphere” at608, i.e. the transmittance and downwelling radiance incident on thetarget, for use in estimating each candidate reflectance at step 610,under the given atmospheric conditions, using surface reflectionparameters for the candidate material obtained from BRDF database 602.The estimated background independent atmospheric correction and theradiance of the image 612 are used to determine an estimate ofbackground reflectance at step 614. The background independent pathterms 606 and the estimated background reflectance from step 614 arethen used at step 616 to determine an estimate of the backgrounddependent atmospheric correction for use in the determination of atleast the candidate spectral signatures. For the reflectance domainprocess shown in FIG. 6, this atmospheric correction is used todetermine a predicted candidate reflectance and error at step 618 and anestimate target reflectance error at step 620 using the radiance fromthe image radiance 612 of the target. For the radiance domain process ofFIG. 7, the estimated background dependent atmospheric correction isused at step 702 to predict the radiance of each candidate materialusing the estimated candidate reflectance from step 610. The targetradiance signature and error is estimated at step 704 from the imageradiance 612. The atmospheric conditions do not need to be taken intoaccount when estimating a target radiance signature.

Referring now to both FIGS. 6 and 7, the estimated target reflectanceand the predicted candidate reflectance for an image constitute,respectively, target and candidate signatures for the reflectancedomain. Similarly, the estimated target radiance and predicted candidateradiance constitute, respectively, target and candidate spectralsignatures for the radiance domain. Target signatures from multipleimages are concatenated, as are candidate signatures for eachacquisition angle, for comparison. Differences between the target andcandidate signatures and error covariance are calculated at step 622. Achi-square test is performed to determine whether the difference issmall enough to declare that the candidate material is a possible matchat step 624. The chi-square test is just one example of a test fordetermining whether there is a match.

Additional details of exemplary methods useful for performing varioussteps of the foregoing processes are given below, and in the sectionsthat follow. Some of these methods, or substantially similar methods,can be adapted, as necessary, for use for certain steps in theimplementation of the material identification process of FIG. 2 forother imaging modalities, which are described below in connection withFIGS. 8-13. All of the methods are carried out by programs executing oncomputing systems. Bold-faced characters indicate vectors.

The aerosol parameters in this example consist of extinctioncoefficients {circumflex over (ε)}={circumflex over (ε)}_(λ), absorptioncoefficients {circumflex over (α)}={circumflex over (α)}_(λ), andasymmetry parameters {circumflex over (ψ)}={circumflex over (ψ)}_(λ),and are first used at step 606 to estimate a set of nominal BIP(background-independent path) terms and three additional setscorresponding to error perturbations δε, δα and δ_(ψ) in {circumflexover (ε)}, {circumflex over (α)} and {circumflex over (ψ)}. The foursets of BIP terms are denoted

IP [4]={

IP,

IP (δε),

IP (δα),

IP (δψ)}. The

IP [4] terms comprise the estimated background independent atmosphericcorrection. These are used with the background radiance L^(b) toestimate the nominal band-integrated background reflectance at step 614and four perturbations denoted {circumflex over (ρ)}[5]={{circumflexover (ρ)}^(b),{circumflex over (ρ)}^(b)(δε),{circumflex over(ρ)}^(b)(δα),{circumflex over (ρ)}^(b)(δψ),{circumflex over(ρ)}^(b)(δL^(b))}. The

IP[4] and {circumflex over (ρ)}^(b) [5] terms are then used at step 616to estimate one set of nominal BDP (background-dependent path) terms andfour additional sets

DP[5]={

DP,

DP (δε),

DP (δα),

DP (δψ),

DP (δL^(b))}.

The

DP [5] terms and the target radiance are then used at step 620 toestimate the band-integrated Lambertian-equivalent target reflectancespectrum in each image and their perturbations, denoted by {circumflexover (ρ)}^(t)[5]={{circumflex over (ρ)}^(t),{circumflex over(β)}^(t)(δε),{circumflex over (ρ)}^(t)(δα),{circumflex over(ρ)}^(t)(δψ), {circumflex over (ρ)}^(t)(δL^(b))}.

The estimated aerosol parameters and background reflectance values{circumflex over (ρ)}^(b) [5] are also used at step 608 to estimate oneset of nominal NEF atmosphere parameters and four error perturbationsdenoted by

TM[5]={

TM,

TM(δε),

TM(δα),

TM(δψ),

TM(δL^(b))}.

The

TM [5] terms are then used at step 610 to estimate the nominal candidatehemispherical directional reflectance and four error perturbations,{circumflex over (ρ)}_(D) ^(c)[5]={{circumflex over (ρ)}_(D) ^(c),{circumflex over (ρ)}_(D) ^(c)(δε), {circumflex over (ρ)}_(D) ^(c)(δα),{circumflex over (ρ)}_(D) ^(c)(δψ), {circumflex over (ρ)}_(D)^(c)(δL^(b))}, and the nominal Lambertian-equivalent solar reflectancesand four error perturbations {circumflex over (ρ)}_(SL)^(c)[5]={{circumflex over (ρ)}_(SL) ^(c), {circumflex over (ρ)}_(SL)^(c)(δε), {circumflex over (ρ)}_(SL) ^(c)(δα), {circumflex over(ρ)}_(SL) ^(c)(δψ), {circumflex over (ρ)}_(SL) ^(c)(δL^(b))}.

The {circumflex over (ρ)}_(D) ^(c)[5], {circumflex over (ρ)}_(SL)^(c)[5] and

DP[5] terms are used to predict at step 618 the band-integratedLambertian-equivalent reflectance spectrum and four error perturbations{circumflex over (ρ)}^(c)[5]={{circumflex over (ρ)}^(c), {circumflexover (ρ)}^(c)(δε), {circumflex over (ρ)}^(c)(δα), {circumflex over(ρ)}^(c)(δψ), {circumflex over (ρ)}^(c)(δL^(b))} of a NEF candidatematerial from BRDF database 602, for each image acquisition geometry.

The {circumflex over (ρ)}^(t)[5] and {circumflex over (ρ)}^(c)[5] termsare used at step 622 to estimate the difference spectrum ξ={circumflexover (ρ)}^(t)−{circumflex over (ρ)}^(c) and four error perturbationsξ[5]={ξ, ξ(δε), ξ(δα), ξ(δψ), ξ(δL^(b))}. The perturbations are used tocalculate the uncertainty in {circumflex over (ρ)}^(t)−{circumflex over(ρ)}^(c) attributable to errors in ε, α, ψ and L^(b). The uncertaintyattributable to errors in target radiance is also calculated, but anerror perturbation is not required because the error does not affect thecandidate reflectance. Consequently, the analytical approach is used toevaluate this term. The uncertainty associated with errors in the NEFmeasurements are also calculated, but the analytical approach is usedbecause these errors do not affect the calculation of the targetreflectance.

The difference spectra are concatenated into a multi-angle signature, anexample of which is shown below, and the error covariance of thissignature is calculated for using in assessing the differences at step622. A chi-square test is performed at step 624 to determine whether thedifferences are small enough to declare that the candidate material is apossible match. Other tests could be used.

$\begin{matrix}{\mspace{79mu}{{\hat{ɛ},\hat{\alpha},\left. \hat{\psi}\rightarrow{\overset{\bullet}{B}{IP}} \right.,{\overset{\bullet}{B}{{IP}({\delta ɛ})}},{\overset{\bullet}{B}{{IP}({\delta\alpha})}},{\overset{\bullet}{B}{{IP}({\delta\psi})}}}\mspace{79mu}{{{\hat{L}}^{b}\overset{\overset{\bullet}{B}{{IP}{\lbrack 4\rbrack}}}{\longrightarrow}\mspace{11mu}{\hat{\rho}}^{b}},{{\hat{\rho}}^{b}({\delta ɛ})},{{\hat{\rho}}^{b}({\delta\alpha})},{{\hat{\rho}}^{b}({\delta\psi})},{{\hat{\rho}}^{b}\left( {\delta\; L^{b}} \right)}}\mspace{79mu}{{\overset{\bullet}{B}{{{IP}\lbrack 4\rbrack}\overset{{\hat{\rho}}^{b}{\lbrack 5\rbrack}}{\longrightarrow}\mspace{14mu}\overset{\bullet}{B}}{DP}},{\overset{\bullet}{B}{{DP}({\delta ɛ})}},{\overset{\bullet}{B}{{DP}({\delta\alpha})}},{\overset{\bullet}{B}{{DP}({\delta\psi})}},{\overset{\bullet}{B}{{DP}\left( {\delta\; L^{b}} \right)}}}{\hat{ɛ},\hat{\alpha},\hat{\psi},\left. {{\hat{\rho}}^{b}\lbrack 5\rbrack}\rightarrow{\overset{\bullet}{A}{TM}} \right.,{\overset{\bullet}{A}{{TM}({\delta ɛ})}},{\overset{\bullet}{A}{{TM}({\delta\alpha})}},{\overset{\bullet}{A}{{TM}({\delta\psi})}},{\overset{\bullet}{A}{{TM}\left( {\delta\; L^{b}} \right)}}}\mspace{79mu}{{{\hat{L}}^{t}\overset{{BDP}{\lbrack 5\rbrack}}{\longrightarrow}\mspace{11mu}{\hat{\rho}}^{t}},{{\hat{\rho}}^{t}({\delta ɛ})},{{\hat{\rho}}^{t}({\delta\alpha})},{{\hat{\rho}}^{t}({\delta\psi})},{{\hat{\rho}}^{t}\left( {\delta\; L^{b}} \right)}}\mspace{79mu}{{\rho_{D}^{c}\lbrack 5\rbrack},{{\rho_{SL}^{c}\lbrack 5\rbrack}\overset{\overset{\bullet}{B}{{DP}{\lbrack 5\rbrack}}}{\longrightarrow}\mspace{11mu}{\hat{\rho}}^{c}},{{\hat{\rho}}^{c}({\delta ɛ})},{{\hat{\rho}}^{c}({\delta\alpha})},{{\hat{\rho}}^{c}({\delta\psi})},{{\hat{\rho}}^{c}\left( {\delta\; L^{b}} \right)}}\mspace{79mu}{{{\hat{\rho}}^{t}\lbrack 5\rbrack},\left. {{\hat{\rho}}^{c}\lbrack 5\rbrack}\rightarrow\xi \right.,{\xi({\delta ɛ})},{\xi({\delta\alpha})},{\xi({\delta\psi})},{\xi\left( {\delta\; L^{b}} \right)}}}} & (4)\end{matrix}$Aerosol Characterization

The aerosol parameters 604 used to estimate the background independentatmospheric correction at step 606 and the BRDF atmosphere at step 608are obtained as follows. The MODTRAN® computer program for modelingatmospheric propagation of electromagnetic radiation contains built inaerosol types. However, in the methods disclosed herein, the effects arepreferably characterized using AERONET (Aerosol Robotic Network) aerosolparameters instead of the aerosol types built into the MODTRAN® program.A cluster analysis study¹ of AERONET data showed that there areeffectively six (6) types of aerosols. The study report listed singlescattering albedo ssa (0.673), extinction coefficient ext (0.673) andasymmetry parameters asym (0.673) at 673 nm. Data at the other threeAERONET wavelengths were obtained from the senior author of the study.These data are listed in Table 2. ¹ Omar, et al., “Development of globalaerosol models using cluster of Aerosol Robotic Network (AERONET)measurements”, JGR, vol. 110, March 2005

TABLE 2 Type 1: Type 2: Type 3: Type 4: Type 5: Type 6: Desert BiomassRural Industrial Polluted Dirty Dust Burning (Background) PollutionMarine Pollution Single scattering albedo @ 441 nm 0.923452 0.8414790.906909 0.937582 0.931900 0.752894 Single scattering albedo @ 673 nm0.927965 0.800305 0.877481 0.921652 0.925049 0.717840 Single scatteringalbedo @ 873 nm 0.926163 0.765592 0.853717 0.905481 0.917595 0.667695Single scattering albedo @ 1022 nm 0.927608 0.745263 0.846585 0.8955030.910512 0.633196 Extinction coefficient @ 441 nm 0.406217 0.3482150.067341 0.370811 0.194477 0.180477 Extinction coefficient @ 673 nm0.326517 0.190405 0.035950 0.190634 0.139910 0.100125 Extinctioncoefficient @ 873 nm 0.291773 0.132512 0.025993 0.125879 0.1145920.071986 Extinction coefficient @ 1022 nm 0.280223 0.111652 0.0228310.101389 0.102874 0.063466 Asymmetry factor @ 441 nm 0.692541 0.6643400.649253 0.694172 0.740076 0.679427 Asymmetry factor @ 673 nm 0.6679370.603481 0.580048 0.611590 0.710529 0.594261 Asymmetry factor @ 873 nm0.672450 0.582214 0.571153 0.574028 0.715632 0.571046 Asymmetry factor @1022 nm 0.673933 0.579726 0.574551 0.562794 0.714759 0.565743

In order to use the AERONET data in the MODTRAN® program for estimatingthe background independent atmospheric correction at step 606 and theBRDF atmosphere at step 608, the AERONET extinction coefficients ext (λ)and single scattering albedos ssa (λ) must first be converted into theMODTRAN normalized extinction coefficients K_(ext) (λ) and normalizedabsorption coefficients K_(abs) (λ). First the absorption coefficientabs (λ) is derived from ext (λ) and ssa (λ) in accordance with the lastline of (1):

$\begin{matrix}{\begin{matrix}{{{ext}(\lambda)} = {{{abs}(\lambda)} + {{scatter}(\lambda)}}} \\{= \left. {{{abs}(\lambda)} + {{{ext}(\lambda)} \cdot {{ssa}(\lambda)}}}\Rightarrow \right.}\end{matrix}{{{abs}(\lambda)} = {{{ext}(\lambda)}\left( {1 - {{ssa}(\lambda)}} \right)}}} & (1)\end{matrix}$

Next, the extinction coefficient at 0.55 microns is estimated by linearinterpolation as in equation (2):

$\begin{matrix}{{{ext}(0.550)} = {{\frac{0.673 - 0.550}{0.673 - 0.441}{{ext}(0.441)}} + {\frac{0.550 - 0.441}{0.673 - 0.441}{{ext}(0.673)}}}} & (2)\end{matrix}$

Finally, K_(ext)(λ) and K_(abs) (λ) are calculated as in equation (3):

$\begin{matrix}{{{K_{ext}(\lambda)} = \frac{{ext}(\lambda)}{{ext}(0.550)}}{{K_{abs}(\lambda)} = \frac{{abs}(\lambda)}{{ext}(0.550)}}} & (3)\end{matrix}$

Since no aerosol vertical profile information is available, and becausethe aerosol effects are only weakly influenced by aerosol verticaldistribution, aerosol is modeled only in the boundary layer.Consequently, the AERONET aerosol parameters are incorporated into theMODTRAN card deck in the following way. For Card2, the ‘ARUSS’ field isset to ‘USS’ to indicate that user-defined aerosol properties will beread from Cards 2D, 2D1 and 2D2. For card2D, IREG(1) is set to 4,indicating that optical properties for the boundary layer aerosols willbe given at 4 wavelengths. IREG(2), IREG(3) and IREG(4) are set to 0indicating that no changes will be made to the other aerosol regimes.For Card2D1, AWCCON is set to blank. For Card2D2, two cards must bewritten, the first providing the aerosol optical properties for thegiven AERONET aerosol type at the first three AERONET wavelengths, andthe second card at the fourth wavelength. Card2D2 #1 specifies, for eachof the first three AERONET wavelengths in Table 6-1, the wavelength λ,in microns, K_(ext)(λ), K_(abs) (λ) and the AERONET asymmetry factorasym(λ). This card therefore has the contents: 0.441, K_(ext)(0.441),K_(abs) (0.441), asym(0.441), 0.673, K_(ext)(0.673), K_(abs) (0.673),asym(0.673), 0.873, K_(ext)(0.873), K_(abs)(0.873), asym(0.873). Thecontents of Card2D2 #2, for the fourth AERONET wavelength in Table 6-1,include 1.022, K_(ext)(1.022), K_(abs) (1.022) and asym(1.022).

From MODTRAN® or other atmospheric modeling program, the followingestimated aerosol spectral parameters {circumflex over (ε)}={circumflexover (ε)}_(λ), {circumflex over (α)}={circumflex over (α)}_(λ) and{circumflex over (ψ)}={circumflex over (ψ)}_(λ) are obtained. It will beassumed for simplicity that errors in these parameters are zero-meanGaussian random variables with known standard deviations that areproportional to the parameter value by factors σ_(δε), σ_(δα) and σ_(δψ)which are the same in each band. It will also be assumed that theseerrors are independent parameter-to-parameter and band-to-band. Thealgorithm can be easily modified to accommodate errors that arecorrelated parameter-to-parameter and band-to-band if such correlationcoefficients are available.

The target radiance signature L^(b) and background radiance signatureL^(t) for each image are assumed to be provided from a multispectralsensor with N_(b) bands. It will be assumed that errors in eachcomponent of these radiance measurements are zero-mean Gaussian randomvariables with standard deviations that are proportional to the radianceby a factor σ_(cal). It will also be assumed that these errors areindependent band-to-band and image-to-image.

Estimation of Spectral BIP Terms

The following is a description of one method for estimating spectral BIPterms.

The spectral BIP terms appropriate for the reflective region consist of,in this example, the following quantities indexed by wavenumber ν. Thesequantities are calculated as functions of aerosol properties andacquisition geometry:

L_(ν) ^(AS)—Solar photons that are scattered into the sensor's field ofview via single or multiple scatter events within the atmosphere withoutever reaching the target or background.

L_(ν) ^(DSR)—Solar photons that pass directly through the atmosphere toa 100% reflectance Lambertian target, reflect off the target, andpropagate directly through the atmosphere to the sensor. Light for thispath term is attenuated by absorption and by light scattering out of thepath, but no radiance is scattered into the path. This path does notinvolve any interaction with the background.

L_(ν) ^(SSR)—Solar photons that are scattered by the atmosphere onto a100% reflectance Lambertian target, reflect off the target, andpropagate directly through the atmosphere to the sensor. This path termdoes not involve any interactions with the background.

L_(ν) ^(BDSR)—Solar photons that are directly transmitted to a 100%reflectance Lambertian background, reflect off the background once, andthen scatter into the field of view of the sensor.

L_(ν) ^(BSSR)—Solar photons that are scattered by the atmosphere atleast once, reflect off a 100% reflectance Lambertian background once,and then scatter into the field of view of the sensor.

S_(ν)—Spherical albedo of the bottom of the atmosphere, which can bethought of as a reflection coefficient of the atmosphere.

The spectral BIP terms can be used to calculate the spectral targetradiance L_(ν) as shown in equation (5), where ρ_(ν) ^(t) is the NEFspectral AEV of the target reflectance and ρ_(ν) ^(b) is the spectralAEV of the background reflectance.

$\begin{matrix}\begin{matrix}{L_{v} = {{\rho_{v}^{SL}L_{v}^{DSR}} +}} & {\ldots\mspace{14mu}{target}\mspace{14mu}{reflected}\mspace{14mu}{solar}\mspace{14mu}{radiance}} \\{{\frac{\rho_{v}^{D}}{1 - {\rho_{v}^{b}S_{v}}}\left( {{\rho_{v}^{b}S_{v}L_{v}^{DSR}} + L_{v}^{SSR}} \right)} +} & {\ldots\mspace{14mu}{target}\mspace{14mu}{scattered}\mspace{14mu}{downwell}} \\{L_{v}^{AS} + {\frac{\rho_{v}^{b}}{1 - {\rho_{v}^{b}S_{v}}}\left( {L_{v}^{BDSR} + L_{v}^{BSSR}} \right)}} & {{{\ldots\mspace{14mu}{atmospheric}}\&}\mspace{11mu}{background}\mspace{14mu}{scatter}}\end{matrix} & (5)\end{matrix}$

Instead of the preceding equation, the following, faster band integratedversion given in equation (6) can be used:L _(i) ^(t)=ρ_(i) ^(SL) L _(i) ^(TRDS)+ρ_(i) ^(D) L _(i) ^(TRH) +L _(i)^(PSMS)  (6)

Three additional sets of spectral BIP terms are calculated, one for eachof the sets of aerosol parameters (δ+Δε,α,ψ), (ε,α+Δα,ψ) and (α, α,ψ+Δψ), where αε, Δα and Δψ are small perturbations. These perturbedterms are used to calculate. Contributions to the uncertainty in theestimated background reflectance attributable to errors in ε, α and ψand contributions to the uncertainty of the band-integrated BDP terms,as explained below.

The ordering of the estimated spectral BIP terms

IP[4] is shown in equation (7):

IP[0]=BIP({circumflex over (ε)},{circumflex over (α)},{circumflex over(ψ)}) . . . nominal BIP terms

IP[1]=({circumflex over (ε)}+Δε,{circumflex over (α)},{circumflex over(ψ)}) . . . extinction perturbation

IP[2]=BIP({circumflex over (ε)},{circumflex over (α)}+Δα,{circumflexover (ψ)}) . . . absorption perturbation

IP[3]=BIP({circumflex over (ε)},{circumflex over (α)},{circumflex over(ψ)}+Δψ) . . . asymmetry perturbation  (7)

A function {circumflex over (ρ)}^(b)(ε, α, ψ, L^(b)) will be implicitlydefined which estimates the band-effective Lambertian-equivalentbackground reflectance {circumflex over (ρ)}_(i) ^(b)={circumflex over(ρ)}^(b)*ε_(i), α_(i), ψ_(i), L_(i) ^(b)) in band i. Since thecalculations are the same in each band, the band index will be droppedin the following equations. First, define the forward model functionL(ρ^(b)) by equation (8), which gives the background radiance for agiven background reflectance ρ^(b) and spectral BIP terms calculatedfrom aerosol parameters ε, α and ψ:

$\begin{matrix}{{L\left( \rho^{b} \right)} = {\int{\left\lbrack {L_{v}^{AS} + {\frac{\rho^{b}}{1 - {\rho^{b}S_{v}}}\left( {L_{v}^{DSR} + L_{v}^{SSR} + L_{v}^{BDSR} + L_{v}^{BSSR}} \right)}} \right\rbrack R_{v}{\mathbb{d}v}}}} & (8)\end{matrix}$

The estimated background reflectance is the value {circumflex over(ρ)}^(b) that solves equation (9).L({circumflex over (ρ)}^(b))={circumflex over (L)} ^(b)  (9)

Equation (9) therefore defines the desired function {circumflex over(ρ)}^(b)(L^(b), ε, α, ψ) by the implicit function theorem. Sinceequation (9) is non-linear in {circumflex over (ρ)}^(b), it is solved byapplying the Newton-Raphson method to equation (10):F({circumflex over (ρ)}^(b))=L({circumflex over (ρ)}^(b))−{circumflexover (L)} ^(b)=0  (10)

The derivative of F required by the Newton-Raphson method is given inequation (11):

$\begin{matrix}\begin{matrix}{{F^{\prime}\left( \rho^{b} \right)} = \frac{\partial L}{\partial\rho^{b}}} \\{= {\int{\frac{1}{\left( {1 - {\rho^{b}S_{v}}} \right)^{2}}\left( {L_{v}^{DSR} + L_{v}^{SSR} + L_{v}^{BDSR} + L_{v}^{BSSR}} \right)R_{v}{\mathbb{d}v}}}}\end{matrix} & (11)\end{matrix}$

The background reflectance and its perturbations {circumflex over(ρ)}^(b)[5] corresponding to perturbations Δε, Δα, Δψ and ΔL^(b) arecalculated for use in calculating the perturbed BDP terms and NEFatmospheres. The notation and ordering are shown in equation (12):{circumflex over (ρ)}^(b)[0]=ρ^(b)(

IP[0],{circumflex over (L)} ^(b)) . . . nominal reflectance{circumflex over (ρ)}^(b)[1]=ρ^(b)(

IP[1],{circumflex over (L)} ^(b)) . . . extinction perturbation{circumflex over (ρ)}^(b)[2]=ρ^(b)(

IP[2],{circumflex over (L)} ^(b)) . . . absorption perturbation{circumflex over (ρ)}^(b)[3]=ρ^(b)(

IP[3],{circumflex over (L)} ^(b)) . . . asymmetry perturbation{circumflex over (ρ)}[4]=ρ^(b)(

IP[0],{circumflex over (L)} ^(b) +ΔL ^(b)) . . . background radianceperturbation  (12)Estimation of Background Reflectance

Following is a description of one example of a method for estimatingbackground reflectance in an image.

The uncertainty contributions to background reflectance corresponding toerrors Δε, Δα, Δψ and ΔL^(b) are calculated as in equation (13):δ{circumflex over (ρ)}^(b)[0]={circumflex over (ρ)}^(b)[1]−{circumflexover (ρ)}^(b)[0] . . . extinction errorδ{circumflex over (ρ)}^(b)[1]={circumflex over (ρ)}^(b)[2]−{circumflexover (ρ)}^(b)[0] . . . absorption errorδ{circumflex over (ρ)}^(b)[2]={circumflex over (ρ)}^(b)[3]−{circumflexover (ρ)}^(b)[0] . . . asymmetry errorδ{circumflex over (ρ)}^(b)[3]={circumflex over (ρ)}^(b)[4]−{circumflexover (ρ)}^(b)[0] . . . background radiance error  (13)

The standard deviation σ_(δ{circumflex over (ρ)}) _(b) _(,δε) of theerror in {circumflex over (ρ)}^(b) attributable to errors in ε isderived in equation (14). It is assumed that the standard deviationσ_(δε) of errors in ε is available from the provider of the aerosolparameters. The last two equations in (14) forσ_(δ{circumflex over (ρ)}) _(b) _(,δα) and σ_(δ{circumflex over (ρ)})_(b) _(,δψ) are derived similarly.

$\begin{matrix}\begin{matrix}{\sigma_{{\delta{\hat{\rho}}^{b}},{\delta ɛ}} = {\frac{\partial{\hat{\rho}}^{b}}{\partial ɛ}\delta_{\delta ɛ}}} & {...\mspace{14mu}{extinction}\mspace{14mu}{contribution}} \\{\approx {\left\lbrack {{{\hat{\rho}}^{b}\left( {ɛ + {\Delta\; ɛ}} \right)} - {{\hat{\rho}}^{b}(ɛ)}} \right\rbrack\frac{\sigma_{\delta ɛ}}{\Delta\; ɛ}}} & \\{= {\delta{{\hat{\rho}}^{b}\lbrack 0\rbrack}\frac{\sigma_{\delta ɛ}}{\Delta ɛ}}} & \\{\sigma_{{\delta{\hat{\rho}}^{b}},{\delta\alpha}} \approx {\left\lbrack {{{\hat{\rho}}^{b}\left( {\alpha + {\Delta\alpha}} \right)} - {{\hat{\rho}}^{b}(\alpha)}} \right\rbrack\frac{\sigma_{\delta\alpha}}{\Delta\;\alpha}}} & {\ldots\mspace{14mu}{absorption}\mspace{14mu}{contribution}} \\{\sigma_{{\delta{\hat{\rho}}^{b}},{\delta\psi}} \approx {\left\lbrack {{{\hat{\rho}}^{b}\left( {\psi + {\Delta\psi}} \right)} - {{\hat{\rho}}^{b}(\psi)}} \right\rbrack\frac{\sigma_{\delta\psi}}{\Delta\;\psi}}} & {\ldots\mspace{14mu}{asymmetry}\mspace{14mu}{contribution}}\end{matrix} & (14)\end{matrix}$

Equation (14) is evaluated as shown in the following equation (15):

$\begin{matrix}{{\sigma_{{\delta{\hat{\rho}}^{b}},{\delta ɛ}} \approx {\delta{{\hat{\rho}}^{b}\lbrack 0\rbrack}\frac{\sigma_{\delta ɛ}}{\Delta ɛ}\mspace{14mu}\ldots\mspace{14mu}{extinction}\mspace{14mu}{contribution}}}\text{}{\sigma_{{\delta{\hat{\rho}}^{b}},{\delta\alpha}} \approx {\delta{{\hat{\rho}}^{b}\lbrack 1\rbrack}\frac{\sigma_{\delta\alpha}}{\Delta\alpha}\mspace{14mu}\ldots\mspace{14mu}{absorption}\mspace{14mu}{contribution}}}\text{}{\sigma_{{\delta{\hat{\rho}}^{b}},{\delta\psi}} \approx {\delta{{\hat{\rho}}^{b}\lbrack 2\rbrack}\frac{\sigma_{\delta\psi}}{\Delta\psi}\mspace{14mu}\ldots\mspace{14mu}{asymmetry}\mspace{14mu}{contribution}}}} & (15)\end{matrix}$

The standard deviation σ_(δ{circumflex over (ρ)}) _(b) _(,δL) _(b) ofthe error in {circumflex over (ρ)}^(b) attributable to errors in L^(b)is derived in equation (16). The second equality follows from thecharacterization of the standard deviation of radiance errors as afraction σ_(cal) of the radiance, and the third equality expresses theanswer in terms of

$\frac{\partial L^{b}}{\partial{\hat{\rho}}^{b}},$which is given by equation (16) evaluated at {circumflex over (ρ)}^(b):

$\begin{matrix}\begin{matrix}{\sigma_{{\delta{\hat{\rho}}^{b}},{\delta\; L^{b}}} = {\frac{\partial{\hat{\rho}}^{b}}{\partial L^{b}}\sigma_{\delta\; L^{b}}}} \\{= {\frac{\partial{\hat{\rho}}^{b}}{\partial L^{b}}\sigma_{{cal}^{L^{b}}}}} \\{= {\frac{L^{b}}{\left( \frac{\partial L^{b}}{\partial{\hat{\rho}}^{b}} \right)}\sigma_{cal}}}\end{matrix} & (16)\end{matrix}$

The standard deviation σ_(δ{circumflex over (ρ)}) _(b) of the errors in{circumflex over (ρ)}^(b) is obtained by combining the standarddeviations of the individual contributors in quadrature as in equation(17), because the errors are assumed to be independent.σ_(δ{circumflex over (ρ)}) _(b) =√{square root over(σ_(δ{circumflex over (ρ)}) _(b) ²+σ_(δ{circumflex over (ρ)}) _(b)_(,δε) ²+σ_(δ{circumflex over (ρ)}) _(b) _(,δα)²+σ_(δ{circumflex over (ρ)}) _(b) _(,δψ) ²)}  (17)Estimate of Band-Integrated BDP Terms

Following is a description of one method for estimating band-integratedBDP terms.

The estimated band-integrated BDP terms

DP are calculated from the estimated spectral BIP terms

IP and the estimated band-effective background reflectance {circumflexover (ρ)}_(i) ^(b) as shown in equation (18).

$\begin{matrix}{{{\hat{L}}_{i}^{TRDS} = {\int{{\hat{L}}_{v}^{DSR}{R_{i}(v)}{\mathbb{d}v}}}}{{\hat{L}}_{i}^{TRH} = {\int{\frac{1}{1 - {{\hat{\rho}}_{i}^{b}S_{v}}}\left( {{{\hat{\rho}}_{i}^{b}{\hat{S}}_{v}{\hat{L}}_{v}^{DSR}} + {\hat{L}}_{v}^{SSR}} \right){R_{i}(v)}{\mathbb{d}v}}}}{{\hat{L}}_{i}^{PSMS} = {\int{\left\lbrack {{\hat{L}}_{v}^{AS} + {\frac{{\hat{\rho}}_{i}^{b}}{1 - {{\hat{\rho}}_{i}^{b}S_{v}}}\left( {{\hat{L}}_{v}^{BDSR} + {\hat{L}}_{v}^{BSSR}} \right)}} \right\rbrack{R_{i}(v)}{\mathbb{d}v}}}}} & (18)\end{matrix}$

Equation (19) will be expressed as

DP=BDP(

IP,{circumflex over (ρ)}^(b)) where the terms {circumflex over(L)}_(DSR), {circumflex over (L)}^(SSR), {circumflex over (L)}^(AS),{circumflex over (L)}^(BDSR), {circumflex over (L)}^(BSSR), Ŝ^(DSR) andR are obtained directly from the BIP structure.

Four additional sets of band-integrated BDP terms are calculated, threecorresponding to the aerosol perturbations Δε, Δα, Δψ and theaccompanying background reflectance perturbations {circumflex over(ρ)}^(b)(ε+Δε, α, ψ, L^(b)), {circumflex over (ρ)}^(b)(ε,α+Δα, ψ,L^(b)), {circumflex over (ρ)}^(b)(ε, α, ψ+Δψ, L^(b)), and onecorresponding to the background radiance perturbation {circumflex over(ρ)}^(b) (ε, α, ψ, L^(b)+ΔL^(b)). The ordering is shown in equation(19):

DP[0]=BDP(

IP[0],{circumflex over (ρ)}^(b)[0]) . . . nominal BDP terms

DP[1]=BDP(

IP[1],{circumflex over (ρ)}^(b)[1]) . . . extinction perturbation

DP[2]=BDP(

IP[2],{circumflex over (ρ)}^(b)[2]) . . . absorption perturbation

DP[3]=BDP(

IP[3],{circumflex over (ρ)}^(b)[3]) . . . asymmetry perturbation

DP[4]=BDP(

IP[0],{circumflex over (ρ)}^(b)[4]) . . . background radianceperturbation  (19)Estimating BRDF Atmosphere

The following method is an example of estimating a “BRDF atmosphere”,which comprises atmospheric transmittances and downwelling radiance. TheBRDF atmosphere is used by the NEFDS to predicted candidate reflectancesignatures. The following example assumes that the target surface ishorizontal.

The estimated NEF atmosphere terms

TM are calculated from the estimated aerosol parameters {circumflex over(ε)}, {circumflex over (α)} and {circumflex over (ψ)}, and from theestimated background reflectance {circumflex over (ρ)}^(b). The processwill be denoted

TM=ATM ({circumflex over (ε)}, {circumflex over (α)}, {circumflex over(ψ)}, {circumflex over (ρ)}^(b)). Four additional sets of NEF atmosphereterms are also calculated, three corresponding to the aerosolperturbations Δε, Δα, Δψ and accompanying background reflectanceperturbations {circumflex over (ρ)}^(b)[1]={circumflex over(ρ)}^(b)(ε+Δε, α, ψ, L^(b)), {circumflex over (ρ)}^(b)[2]={circumflexover (ρ)}^(b)(ε, α+Δα, ψ, L^(b)), {circumflex over (ρ)}^(b)(ε, α, ψ+Δψ,L^(b)) and one corresponding to a perturbation of the backgroundradiance {circumflex over (ρ)}^(b)[4]={circumflex over (ρ)}^(b)(ε, α, ψ,ΔL^(b)). These perturbed terms are used to calculate derivatives of theNEF directional hemispherical reflectance and Lambertian-equivalentsolar reflectance AEVs. The ordering is shown in equation (20):

TM[0]=ATM({circumflex over (ε)},{circumflex over (α)},{circumflex over(ψ)},{circumflex over (ρ)}^(b)[0]) . . . nominal BDP terms

TM[1]=ATM({circumflex over (ε)}+Δε,{circumflex over (α)}ε,{circumflexover (ψ)},{circumflex over (ρ)}^(b)[1]) . . . extinction perturbation

TM[2]=ATM({circumflex over (ε)},{circumflex over (α)}+Δ{circumflex over(α)},{circumflex over (ψ)},{circumflex over (ρ)}^(b)[2]) . . .absorption perturbation

TM[3]=ATM({circumflex over (ε)},{circumflex over (α)}+Δ{circumflex over(α)},{circumflex over (ψ)},{circumflex over (ρ)}^(b)[2]) . . . asymmetryperturbation

TM[4]=ATM({circumflex over (ε)},{circumflex over (α)},{circumflex over(ψ)},{circumflex over (ρ)}^(b)[4]) . . . background radianceperturbation  (20)Estimation of Target Reflectance

The following method is an example for estimating target reflectancewith the process of FIG. 6, as well as material identification processesfor different imaging modalities that require estimation of targetreflectance from a radiance measurement. This example assumes that theorientation of the target surface is horizontal.

The band-effective Lambertian-equivalent target reflectance ρ_(i) ^(t)in band i is estimated using equation (21), where the second linefollows from the Lambertian-equivalent assumption ρ_(i) ^(SL)=ρ_(i)^(D)=ρ_(i) ^(t):

$\begin{matrix}\begin{matrix}{L_{i}^{t} = {{\rho_{i}^{SL}L_{i}^{TRDS}} + {\rho_{i}^{D}L_{i}^{TRH}} + L_{i}^{PSMS}}} \\{= {{\rho_{i}^{t}\left( {L_{i}^{TRDS} + L_{i}^{TRH}} \right)} + L_{i}^{PSMS}}}\end{matrix} & (21)\end{matrix}$

Equation (21) is now solved for ρ_(i) ^(t) as in equation (22), where{circumflex over (L)}^(s)={circumflex over (L)}^(PSMS) is the radiancescattered by the atmosphere and background, and {circumflex over(L)}^(r)={circumflex over (L)}^(PSMS) radiance reflected by a perfectLambertian reflector.

$\begin{matrix}{{\hat{\rho}}^{t} = \frac{{\hat{L}}^{t} - {\hat{L}}^{s}}{{\hat{L}}^{r}}} & (22)\end{matrix}$

Equation (22) will be expressed as {circumflex over (ρ)}=ρ^(t)(

DP), where {circumflex over (L)}^(s) and {circumflex over (L)}^(r) arecalculated from the

DP terms {circumflex over (L)}^(PSMS), {circumflex over (L)}^(TRDS) and{circumflex over (L)}^(TRH) explained above.

Four additional reflectance perturbations are also calculated as inequation (23), corresponding to the three perturbations of theband-integrated BDP terms and the perturbation in the backgroundradiance. These perturbations will be used to calculate the correlatederror in the delta reflectance and uncertainty as described below.{circumflex over (ρ)}^(t)[0]=ρ^(t)(

DP[0]) . . . nominal reflectance{circumflex over (ρ)}^(t)[1]=ρ^(t)(

DP[1]) . . . extinction perturbation{circumflex over (ρ)}^(t)[2]=ρ^(t)(

DP[2]) . . . absorption perturbation{circumflex over (ρ)}^(t)[3]=ρ^(t)(

DP[3]) . . . asymmetry perturbation{circumflex over (ρ)}^(t)[4]=ρ^(t)(

DP[4]) . . . background radiance perturbation  (23)

The uncertainty contributions to target reflectance corresponding toerrors Δε, Δα, Δψ, as ΔL^(b) are calculated as in equation (24):δ{circumflex over (ρ)}^(t)[0]={circumflex over (ρ)}^(t)[1]−{circumflexover (ρ)}^(t)[0] . . . extinction errorδ{circumflex over (ρ)}^(t)[1]={circumflex over (ρ)}^(t)[2]−{circumflexover (ρ)}^(t)[0] . . . absorption error (24)δ{circumflex over (ρ)}^(t)[2]={circumflex over (ρ)}^(t)[3]−{circumflexover (ρ)}^(t)[0] . . . asymmetry errorδ{circumflex over (ρ)}^(t)[3]={circumflex over (ρ)}^(t)[4]−{circumflexover (ρ)}^(t)[0] . . . background radiance error  (24)

The uncertainty contribution from target radiance errors is calculatedas in equation (25). The first equality follows from first-order errorpropagation. The first term in the second equality is obtained bydifferentiating equation (22), and the second term follows from thecharacterization of radiance error as a percentage σ_(cal) of theradiance:

$\begin{matrix}\begin{matrix}{\sigma_{{\delta{\hat{\rho}}^{t}},{\delta\; L^{t}}} = {{\frac{\partial{\hat{\rho}}^{t}}{\partial L^{t}}}\sigma_{\delta\; L^{t}}}} \\{= {{{- \frac{1}{L^{r}}}}\sigma_{cal}L^{t}}} \\{= {\sigma_{cal}\frac{L^{t}}{L^{r}}}}\end{matrix} & (25)\end{matrix}$Estimation of Candidate Reflectance

For the process described herein, such as those of FIGS. 6 and 7, theband-effective Lambertian-equivalent reflectance for a candidatematerial can be calculated as shown in equation (26), where {circumflexover (ρ)}^(D) is the estimated AEV directional hemispherical reflectanceand {circumflex over (ρ)}^(SL) is the estimated AEVLambertian-equivalent solar reflectance calculated by the NEF using theestimated NEF atmosphere.

$\begin{matrix}\begin{matrix}{{\hat{\rho}}^{c} = {{\frac{{\hat{L}}^{TRH}}{{\hat{L}}^{TRH} + {\hat{L}}^{TRDS}}{\hat{\rho}}^{D}} + {\frac{{\hat{L}}^{TRDS}}{{\hat{L}}^{TRH} + {\hat{L}}^{TRDS}}{\hat{\rho}}^{SL}}}} \\{= {{\hat{f} \cdot {\hat{\rho}}^{D}} + {\left( {1 - \hat{f}} \right) \cdot {\hat{\rho}}^{SL}}}}\end{matrix} & (26)\end{matrix}$

Equation (26) will be expressed as {circumflex over (ρ)}^(c)=ρ^(c)(

DP, {circumflex over (ρ)}^(D), {circumflex over (ρ)}^(SL)), where{circumflex over (f)} is calculated from the BDP terms L^(TRH) andL^(TRDS) as explained above.

Four additional reflectance perturbations are also calculated as inequation (27), corresponding to the three perturbations of theband-integrated BDP terms and a perturbation in the backgroundreflectance. These perturbation will be used to calculate the correlatederror in the delta reflectance.{circumflex over (ρ)}^(c)[0]=ρ^(c)({circumflex over (ρ)}_(D)^(c)[0],{circumflex over (ρ)}_(SL) ^(c)[0] . . . nominal reflectance{circumflex over (ρ)}^(c)[1]=ρ^(c)(

DP[1],{circumflex over (ρ)}_(D) ^(c)[1],{circumflex over (ρ)}_(SL)^(c)[1] . . . extinction perturbation{circumflex over (ρ)}^(c)[2]=ρ^(c)(

DP[2],{circumflex over (ρ)}_(D) ^(c)[2],{circumflex over (ρ)}_(SL)^(c)[2] . . . absorption perturbation{circumflex over (ρ)}^(c)[3]=ρ^(c)(

DP[3],{circumflex over (ρ)}_(D) ^(c)[3],{circumflex over (ρ)}_(SL)^(c)[3]) . . . asymmetry perturbation{circumflex over (ρ)}^(c)[4]=ρ^(c)(

DP[4],{circumflex over (ρ)}_(D) ^(c)[4],{circumflex over (ρ)}_(SL)^(c)[4] . . . background radiance perturbation  (27)

The uncertainty contributions to candidate reflectance corresponding toerrors Δε, Δα, Δψ and ΔL^(b) are calculated as in equation (28):δ{circumflex over (ρ)}^(c)[0]={circumflex over (ρ)}^(c)[1]−{circumflexover (ρ)}^(c)[0] . . . extinction errorδ{circumflex over (ρ)}^(c)[1]={circumflex over (ρ)}^(c)[2]−{circumflexover (ρ)}^(c)[0] . . . absorption errorδ{circumflex over (ρ)}^(c)[2]={circumflex over (ρ)}^(c)[3]−{circumflexover (ρ)}^(c)[0] . . . asymmetry errorδ{circumflex over (ρ)}^(c)[3]={circumflex over (ρ)}^(c)[4]−{circumflexover (ρ)}^(c)[0] . . . background radiance error  (28)

An uncertainty contribution attributable to NEF BRDF measurement errorsis calculated in accordance with equation (29):

$\begin{matrix}\begin{matrix}{\sigma_{{\delta\rho}^{c},{\delta\;\rho^{NEF}}}^{2} = {{\begin{bmatrix}\frac{\partial{\hat{\rho}}^{c}}{\partial{\hat{\rho}}^{D}} & \frac{\partial{\hat{\rho}}^{c}}{\partial{\hat{\rho}}^{SL}}\end{bmatrix}\begin{bmatrix}\sigma_{\rho^{D}}^{2} & \sigma_{\rho^{D},\rho^{SL}} \\\sigma_{\rho^{SL},\rho^{D}} & \sigma_{\rho^{SL}}^{2}\end{bmatrix}}\begin{bmatrix}\frac{\partial{\hat{\rho}}^{c}}{\partial{\hat{\rho}}^{D}} \\\frac{\partial{\hat{\rho}}^{c}}{\partial{\hat{\rho}}^{SL}}\end{bmatrix}}} \\{= {{\begin{bmatrix}\hat{f} & {1 - \hat{f}}\end{bmatrix}\begin{bmatrix}\sigma_{\rho^{D}}^{2} & \sigma_{\rho^{D},\rho^{SL}} \\\sigma_{\rho^{SL},\rho^{D}} & \sigma_{\rho^{SL}}^{2}\end{bmatrix}}\begin{bmatrix}\hat{f} \\{1 - \hat{f}}\end{bmatrix}}} \\{= {{{\hat{f}}^{2}\sigma_{\rho^{D}}^{2}} + {2{\hat{f}\left( {1 - \hat{f}} \right)}\sigma_{\rho^{D},\rho^{SL}}} + {\left( {1 - \hat{f}} \right)^{2}\sigma_{\rho^{SL}}^{2}}}}\end{matrix} & (29)\end{matrix}$Delta Reflectance and Uncertainty

One method of performing step 622 of FIG. 6 is, as mentioned above,determining delta reflectance and uncertainty. One method for makingthis determination is described below. Similar methods can be used inconnection with implementations for other imaging modalities describedbelow.

Turning first to a step 622, for a single image case, start withX_(i)=(ε_(i), α_(i), ψ_(t), L_(i) ^(b)), X=(X₁, . . . , X_(N) _(b) ),Y_(i)=(ρ_(i) ^(D),ρ_(i) ^(SL)) and Y=(Y₁, . . . , Y_(N) _(b) . The deltareflectance spectrum is given by equation (30):

$\begin{matrix}\begin{matrix}{\xi = {\xi\left( {X,L^{t},Y} \right)}} \\{= {{{\hat{\rho}}^{t}\left( {X,L^{t}} \right)} - {{\hat{\rho}}^{c}\left( {X,Y} \right)}}}\end{matrix} & (30)\end{matrix}$

Equation (30) implies that errors in {circumflex over (ρ)}^(t) and{circumflex over (ρ)}^(c) will be correlated through a common dependenceon X. In the single image case, errors in X, L^(t), and Y areindependent. Writing the combined variables as Z=(X, L^(t), Y), itfollows that S_(δξ)=cov(δξ) is given by equation (31).

$\begin{matrix}\begin{matrix}{S_{\delta\xi} = {\frac{\partial\xi}{\partial Z}{{cov}\left( {\delta\; Z} \right)}\frac{\partial\xi^{T}}{\partial Z}}} \\{= {{\frac{\partial\xi}{\partial Z}\begin{bmatrix}{{cov}\left( {\delta\; X} \right)} & \; & \; \\\; & {{cov}\left( {\delta\; L^{t}} \right)} & \; \\\; & \; & {{cov}\left( {\delta\; Y} \right)}\end{bmatrix}}\frac{\partial\xi^{T}}{\partial Z}}} \\{= {{\frac{\partial\xi}{\partial X}{{cov}\left( {\delta\; X} \right)}\frac{\partial\xi^{T}}{\partial X}} + {\frac{\partial\xi}{\partial L^{t}}{{cov}\left( {\delta\; L^{t}} \right)}\frac{\partial\xi^{T}}{\partial L^{t}}} + {\frac{\partial\xi}{\partial Y}{{cov}\left( {\delta\; Y} \right)}\frac{\partial\xi^{T}}{\partial Y}}}} \\{= {{\frac{\partial\xi}{\partial X}{{cov}\left( {\delta\; X} \right)}\frac{\partial\xi^{T}}{\partial X}} + {\frac{\partial{\hat{\rho}}^{t}}{\partial L^{t}}{{cov}\left( {\delta\; L^{t}} \right)}\frac{\partial{\hat{\rho}}^{t^{T}}}{\partial L^{t}}} + {\frac{\partial{\hat{\rho}}^{c}}{\partial Y}{{cov}\left( {\delta\; Y} \right)}\frac{\partial{\hat{\rho}}^{c^{T}}}{\partial Y}}}}\end{matrix} & (31)\end{matrix}$

Equation (31) simplifies to the form shown in equation (32). The term

$\frac{\partial\xi}{\partial X}$term is an N_(b)×N_(b) block diagonal matrix, one block per band, with1×4 sub-blocks, because each ξ_(i) depends only on ε_(i), α_(i), ψ_(i)and L_(i) ^(b) and not on ε_(j), α_(j), ψ_(j) or L_(j) ^(b) for j≠i. Thecov(δX) term is an N_(b)×N_(b) block diagonal matrix with 4×4 sub-blocksbecause it is assumed that errors in X are independent band-to-band andparameter-to-parameter. Similarly,

$\frac{\partial{\hat{\rho}}^{t}}{\partial L^{t}}$and cov(δL^(t)) are N_(b)×N_(b) diagonal matrices. The

$\frac{\partial{\hat{\rho}}^{c}}{\partial Y}$term is an N_(b)×N_(b) block diagonal matrix, one block per band, with1×2 sub-blocks, because each {circumflex over (ρ)}_(i) ^(c) depends onlyon {circumflex over (ρ)}_(i) ^(D) and {circumflex over (ρ)}_(i) ^(SL),and not on {circumflex over (ρ)}_(j) ^(D) or {circumflex over (ρ)}_(j)^(SL) for j≠i. The cov(δY) term is an N_(b)×N_(b) block diagonal matrixwith one block per band (because for individual materials the NEF doesnot compute the error covariance between bands), where each sub-block isa 2×2 matrix supplied by the NEF. Altogether, equation (32) shows thatS_(δξ) is an N_(b)×N_(b) diagonal matrix.

$\begin{matrix}\begin{matrix}{S_{\delta\xi} = {{\frac{\partial\xi}{\underset{\underset{\underset{{1 \times 4\mspace{14mu}{sub}} - {blocks}}{N_{b} \times N_{b}{blk}\mspace{14mu}{diag}}}{︸}}{\partial X}}\underset{\underset{\underset{{4 \times 4\mspace{14mu}{sub}} - {blocks}}{N_{b} \times N_{b}{blk}\mspace{14mu}{diag}}}{︸}}{{cov}\left( {\delta\; X} \right)}\frac{\partial\xi^{T}}{\underset{\underset{\underset{{4 \times 1\mspace{14mu}{sub}} - {blocks}}{N_{b} \times N_{b}{blk}\mspace{14mu}{diag}}}{︸}}{\partial X}}\mspace{14mu}\ldots\mspace{14mu} 4 \times 4\mspace{14mu}{diagonal}} +}} \\{{\frac{\partial{\hat{\rho}}^{t}}{\underset{\underset{N_{b} \times N_{b}\mspace{14mu}{diagonal}}{︸}}{\partial L^{t}}}\underset{\underset{N_{b} \times N_{b}\mspace{14mu}{diagonal}}{︸}}{{cov}\left( {\delta\; L^{t}} \right)}\frac{\partial{\hat{\rho}}^{t^{T}}}{\underset{\underset{N_{b} \times N_{b}\mspace{14mu}{diagonal}}{︸}}{\partial L^{t}}}\mspace{14mu}\ldots\mspace{14mu} 4 \times 4\mspace{14mu}{diagonal}} +} \\{{\frac{\partial{\hat{\rho}}^{c}}{\underset{\underset{\underset{{1 \times 2\mspace{14mu}{sub}} - {blocks}}{N_{b} \times N_{b}{blk}\mspace{14mu}{diag}}}{︸}}{\partial Y}}\underset{\underset{\underset{{2 \times 2\mspace{14mu}{sub}} - {blocks}}{N_{b} \times N_{b}{blk}\mspace{14mu}{diag}}}{︸}}{{cov}\left( {\delta\; Y} \right)}\frac{\partial{\hat{\rho}}^{c^{T}}}{\underset{\underset{\underset{{2 \times 1\mspace{14mu}{sub}} - {blocks}}{N_{b} \times N_{b}{blk}\mspace{14mu}{diag}}}{︸}}{\partial Y}}\mspace{14mu}\ldots\mspace{14mu} 4 \times 4\mspace{14mu}{diagonal}} +}\end{matrix} & (32)\end{matrix}$

In evaluating equation (32), the only term that has not already beencalculated is

$\frac{\partial\xi}{\partial X}.$This term has the N_(b)×N_(b) block diagonal form with 1×4 sub-blocksshown in equation (33):

$\begin{matrix}{{\frac{\partial\xi}{\partial X} = \begin{bmatrix}\frac{\partial\xi_{1}}{\partial X_{1}} & \; & \; \\\; & \ddots & \; \\\; & \; & \frac{\partial\xi_{N_{b}}}{\partial X_{N_{b}}}\end{bmatrix}}{\frac{\partial\xi_{i}}{\partial X_{i}} = \left\lbrack {\frac{\partial\xi_{i}}{\partial ɛ_{i}}\frac{\partial\xi_{i}}{\partial\alpha_{i}}\frac{\partial\xi_{i}}{\partial\psi_{i}}\frac{\partial\xi_{i}}{\partial L_{i}^{b}}} \right\rbrack}} & (33)\end{matrix}$

Each entry

$\frac{\partial\xi_{i}}{\partial ɛ_{i}}$in equation (33) is calculated from quantities described in previoussections as shown in equation (34), where the band index i has beendropped:

$\begin{matrix}\begin{matrix}{{\frac{\partial\xi}{\partial ɛ}{\Delta ɛ}} \approx {{\xi\left( {ɛ + {\Delta ɛ}} \right)} - {\xi(ɛ)}}} \\{= {\left\lbrack {{{\hat{\rho}}^{t}\left( {ɛ + {\Delta ɛ}} \right)} - {{\hat{\rho}}^{c}\left( {ɛ + {\Delta ɛ}} \right)}} \right\rbrack - \left\lbrack {{{\hat{\rho}}^{t}(ɛ)} - {{\hat{\rho}}^{c}(ɛ)}} \right\rbrack}} \\{= {\left\lbrack {{{\hat{\rho}}^{t}\left( {ɛ + {\Delta ɛ}} \right)} - {{\hat{\rho}}^{t}(ɛ)}} \right\rbrack - \left\lbrack {{{\hat{\rho}}^{c}\left( {ɛ + {\Delta ɛ}} \right)} - {{\hat{\rho}}^{c}(ɛ)}} \right\rbrack}} \\{= {{\delta{{\hat{\rho}}^{t}({\Delta ɛ})}} - {\delta{{\hat{\rho}}^{c}({\Delta ɛ})}}}} \\{= {{\delta{{\hat{\rho}}^{t}\lbrack 0\rbrack}} - {\delta{{\hat{\rho}}^{c}\lbrack 0\rbrack}}}}\end{matrix} & (34)\end{matrix}$

By the same argument, the other entries of each of the blocks inequation (33) can be calculated as shown in equation (35).

$\begin{matrix}{\mspace{79mu}{{{{\delta\xi}\lbrack 0\rbrack} = {{\frac{\partial\xi}{\partial ɛ}{\Delta ɛ}} \approx {{\delta{{\hat{\rho}}^{t}\lbrack 0\rbrack}} - {\delta{{\hat{\rho}}^{c}\lbrack 0\rbrack}\mspace{14mu}\ldots\mspace{14mu}{extinction}\mspace{14mu}{contribution}}}}}\mspace{79mu}{{{\delta\xi}\lbrack 1\rbrack} = {{\frac{\partial\xi}{\partial\alpha}{\Delta\alpha}} \approx {{\delta{{\hat{\rho}}^{t}\lbrack 1\rbrack}} - {\delta{{\hat{\rho}}^{c}\lbrack 1\rbrack}\mspace{14mu}\ldots\mspace{14mu}{absorption}\mspace{14mu}{contribution}}}}}\mspace{79mu}{{{\delta\xi}\lbrack 2\rbrack} = {{\frac{\partial\xi}{\partial\psi}{\Delta\psi}} \approx {{\delta{{\hat{\rho}}^{t}\lbrack 2\rbrack}} - {\delta{{\hat{\rho}}^{c}\lbrack 2\rbrack}\mspace{14mu}\ldots\mspace{14mu}{asymmetry}\mspace{14mu}{contribution}}}}}{{{\delta\xi}\lbrack 3\rbrack} = {{\frac{\partial\xi}{\partial L^{b}}\Delta\; L^{b}} \approx {{\delta{{\hat{\rho}}^{t}\lbrack 3\rbrack}} - {\delta{{\hat{\rho}}^{c}\lbrack 3\rbrack}\mspace{14mu}\ldots\mspace{14mu}{background}\mspace{14mu}{radiance}\mspace{14mu}{contribution}}}}}}} & (35)\end{matrix}$

Entry i of the diagonal matrix cov (o) is therefore given by equation(36):

$\begin{matrix}{\sigma_{{\delta\xi}_{i}}^{2} = {{\left( \frac{\partial\xi_{i}}{\partial ɛ_{i}} \right)^{2}\sigma_{{\delta ɛ}_{i}}^{2}} + {\left( \frac{\partial\xi_{i}}{\partial\alpha_{i}} \right)^{2}\sigma_{{\delta\alpha}_{i}}^{2}} + {\left( \frac{\partial\xi_{i}}{\partial\psi_{i}} \right)^{2}\sigma_{{\delta\psi}_{i}}^{2}} + {\left( \frac{\partial\xi_{i}}{\partial L_{i}^{b}} \right)^{2}\sigma_{\delta\; L_{i}^{b}}^{2}} + {\left( \frac{\partial{\hat{\rho}}_{i}^{t}}{\partial L_{i}^{t}} \right)^{2}\sigma_{\delta\; L_{i}^{t}}^{2}} + \sigma_{{\delta\;\rho_{i}^{c}},{\delta\rho}_{i}^{NEF}}^{2}}} & (36)\end{matrix}$

Equation (36) is evaluated using equations (25), (29) and (35) as shownin equation (37):

$\begin{matrix}\begin{matrix}{\sigma_{{\delta\xi}_{i}}^{2} = {{\left( {{{\delta\xi}\lbrack 0\rbrack}\frac{\sigma_{\delta ɛ}}{\Delta ɛ}} \right)^{2}\ldots\mspace{14mu}{extinction}\mspace{14mu}{contribution}} +}} \\{{\left( {{{\delta\xi}\lbrack 1\rbrack}\frac{\sigma_{\delta\alpha}}{\Delta\alpha}} \right)^{2}\ldots\mspace{14mu}{absorption}\mspace{14mu}{contribution}} +} \\{{\left( {{{\delta\xi}\lbrack 2\rbrack}\frac{\sigma_{\delta\psi}}{\Delta\psi}} \right)^{2}\ldots\mspace{14mu}{asymmetry}\mspace{14mu}{contribution}} +} \\{{\left( {{{\delta\xi}\lbrack 3\rbrack}\frac{\sigma_{\delta\; L^{b}}}{{\Delta L}^{b}}} \right)^{2}\ldots\mspace{14mu}{background}\mspace{14mu}{radiance}\mspace{14mu}{contribution}} +} \\{{\left( \frac{L^{t}}{L^{r}} \right)\sigma_{cal}^{2}\mspace{14mu}\ldots\mspace{14mu}{target}\mspace{14mu}{radiance}\mspace{14mu}{contribution}} +} \\{\sigma_{{\delta\rho}_{i}^{c},{\delta\rho}_{i}^{NEF}}^{2}\mspace{14mu}\ldots\mspace{14mu}{NEF}\mspace{14mu}{error}\mspace{14mu}{contribution}}\end{matrix} & (37)\end{matrix}$

Turning now to a multiple image case for determining delta reflectanceand uncertainty, start with ξ=[ξ₁ . . . ξ_(N) _(images) ]^(T), whereξ_(j)(Â, {circumflex over (L)}^(b), {circumflex over (L)}^(t),{circumflex over (ρ)}_(j) ^(D), {circumflex over (ρ)}_(j) ^(SL)) is thedelta-reflectance spectrum for image j and Â=[{circumflex over (ε)},{circumflex over (α)}, {circumflex over (ψ)}]^(T) are the estimatedaerosol parameters. Errors in the components of ξ will be correlatedthrough a common dependence on Â, but not on any of the other parametersbecause errors in the measurements for image j (i.e. {circumflex over(L)}^(b), {circumflex over (L)}_(j) ^(T), {circumflex over (ρ)}_(j) ^(D)and {circumflex over (ρ)}_(j) ^(SL)) are uncorrelated with each otherand with those for image k (i.e. {circumflex over (L)}_(k) ^(b),{circumflex over (L)}_(k) ^(t), {circumflex over (ρ)}_(k) ^(D) and{circumflex over (ρ)}_(k) ^(SL)) for j≠k.

The covariance S_(δξ) of errors in ξ therefore has the structure shownin equation (38):

$\begin{matrix}{S_{\delta\xi} = \begin{bmatrix}\ddots & \; & \; & \; & \; \\\; & {{cov}\left( {{\delta\xi}_{j},{\delta\xi}_{j}} \right)} & \; & {{cov}\left( {{\delta\xi}_{j},{\delta\xi}_{k}} \right)} & \; \\\; & \; & \ddots & \; & \; \\\; & {{cov}\left( {{\delta\xi}_{k},{\delta\xi}_{j}} \right)} & \; & {{cov}\left( {{\delta\xi}_{k},{\delta\xi}_{k}} \right)} & \; \\\; & \; & \; & \; & \ddots\end{bmatrix}} & (38)\end{matrix}$

In equation (38) j and k are image numbers. Each cov (δξ_(j), δξ_(j)) isa diagonal matrix whose i^(th) entry is given by equation (36). Each coy(δξ_(j), δξ_(k)) is a diagonal matrix whose i^(th) entry is given byequation (39):

$\begin{matrix}{\sigma_{{\delta\xi}_{j,i},{\delta\xi}_{k,i}} = {{\frac{\partial\xi_{j,i}}{\partial ɛ_{i}}\frac{\partial\xi_{k,i}}{\partial ɛ_{i}}\sigma_{{\delta ɛ}_{i}}^{2}} + {\frac{\partial\xi_{j,i}}{\partial\alpha_{i}}\frac{\partial\xi_{k,i}}{\partial\alpha_{i}}\sigma_{{\delta\alpha}_{i}}^{2}} + {\frac{\partial\xi_{j,i}}{\partial\psi_{i}}\frac{\partial\xi_{k,i}}{\partial\psi_{i}}\sigma_{{\delta\psi}_{i}}^{2}}}} & (39)\end{matrix}$Calculation of Chi-Square Statistic

A chi-square statistic can be determined for the processes described inFIGS. 6 and 7 according to the following method. The method can beadapted for use in the processes described below and for other imagingmodalities. The chi-square statistic is calculated as shown in equation(40) in both the single and multiple image cases. In the single imagecase, ξ is an [N_(bands)×1 vector and S_(δξ) is an N_(bands)·N_(bands)]diagonal matrix whose i^(th) entry is given by equation (36). In themultiple image case ξ is an [N_(images)·N_(bands)]×1 vector and S_(δξ)is an [N_(images)·N_(bands)]×[N_(images)·N_(bands)] matrix given byequations (38), (36) and (39).χ²=ξ^(T) ·S _(δξ) ⁻¹·ξ  (40)Performance of the Chi-Square Test for a Given Probability of Detection

When the correct candidate material is chosen, the expected value of{circumflex over (ρ)}^(t) and {circumflex over (ρ)}^(c) will be equaland the expected value of ξ is 0. In this case, the χ² statistic shouldhave a chi-square distribution with number of degrees of freedomν=N_(bands)·N_(images). Denote its inverse cumulative distributionfunction by Q_(χ) ₂ . For a given probability of detection P_(d), athreshold x is set as shown in equation (41):x=Q _(χ) ₂ (P _(d),ν)  (41)

The chi-square test now reduces to the threshold test shown in equation(42):χ² ≦x

candidate material matches the targetχ² >x

candidate material does not match the target  (42)

This method can be adapted for in connection with the implementationsfor other imaging modalities described in FIGS. 8-13.

Non-Polarimetric Emissive Processing

Non-polarimetric emissive imagery—imagery generated by sensing theemitted and reflected radiance from the physical target using anon-polarimetric sensor—may be a single band; however, likenon-polarimetric reflective imagery, non-polarimetric emissive imageryis preferably either multi-band or hyperspectral. Error sources fornon-polarimetric imagery may include atmospheric conditions, radiometriccalibration, BRDF measurements, and spectral calibration.

FIG. 8 illustrates an example of an implementation of the materialidentification process of FIG. 2 for non-polarimetric emissive imagery.For each image data set, the example in FIG. 8 utilizes the followinginputs: information about the atmosphere through which emitted radiationwas sensed; target and background radiance—such as that used inaccordance with the non-polarimetric reflective case described above;the relative spectral response functions for each band in the imagesdata set; candidate material surface properties; and uncertainties.

Data regarding the atmosphere includes the same aerosol parameters asdescribed above in the non-polarimetric reflective case as well as anair pressure profile, air temperature profile, and humidity profile as afunction of altitude above the target. Candidate material surfaceproperties include directional emissivity and, as previously mentioned,diffuse reflectance and Lambertian-equivalent solar reflectance.Uncertainty inputs include radiance measurement uncertainty, errorcovariance matrix for the diffuse reflectance, Lambertian-equivalentsolar reflectance and directional emissivity, aerosol parameteruncertainty, atmospheric uncertainty, and spectral calibrationuncertainty for hyperspectral sensor imagery.

At step 802, a background-independent atmospheric correction, comprisedof background-independent path terms, is estimated using the aerosolparameters, air pressure, air temperature, and humidity profilesdetermined in accordance with the altitude at which the sensor waslocated above the target when sensing the emissivity. Thebackground-independent atmospheric path terms, shown below in Table 3,are combined with the image radiance data to estimate backgroundreflectance and emissivity of the target as shown at step 804.

TABLE 3 BIP Term Description L_(AE) Atmospherically emitted photons thatpropagate to the sensor aperture L_(AER) Atmospherically emitteddownwelling radiation that is reflected by a 100% Lambertian target andpropagates directly to the sensor aperture L_(BAER) Atmosphericallyemitted downwelling radiation that is reflected from a 100% Lambertianbackground once and scattered into the line of sight T_(AD) Directtransmission of the atmosphere between the target and the sensor T_(AS)Fraction of the background-leaving radiance that is transmitted to thesensor

At step 806, the input data, background reflectance and emissivity, andbackground-independent path terms are combined to estimate abackground-dependent atmospheric correction comprised ofbackground-dependent path terms, shown below in Table 4. At step 808,the image radiance data is used to estimate the target radiance and theerror corresponding to the target radiance.

TABLE 4 BIP Term Description L_(PT) Photons emitted by the atmosphere orbackground that eventually scatter into the line of sight and propagateinto the sensor field of view F_(MS BRDS) Fraction of the backgroundreflected solar radiation that is multiply scattered

At step 810, the process uses the input data to estimate the BRDFatmosphere, which comprises atmospheric transmittances and downwellingradiance on the target. These terms are used to estimate the candidatereflectances and emissivity. At step 812, one or more candidatematerials are selected from BRDF database 813; then, for each candidatematerial, band-effective Lambertian-equivalent candidate reflectance andemissivity spectra are predicted.

At step 814, the best-fitting candidate spectral radiance signature andits error are estimated using the estimated background-dependentatmospheric correction, target radiance signature and error, and thecandidate reflectance and emissivity spectra determined at steps 806,808 and 812, respectively. There are at least two methods forcalculating the best-fitting candidate spectral radiance signature, andeach method is dependent upon the temperature of the multi-angle targetradiance signature calculated in step 808 and a candidate spectralradiance signature determined for one or more particular bands in theimages data set.

FIG. 9 illustrates a first method for calculating the best-fittingcandidate spectral radiance signature, wherein the best-fittingcandidate spectral radiance signature selected is determined byselecting a candidate spectral radiance signature having a temperaturethat best fits the temperature of the multi-angle target radiancesignature. As illustrated in FIG. 9, the ground-level air temperature isprovided as an initial temperature T in step 902. In step 904, acandidate spectral radiance signature and its uncertainty are calculatedusing the temperature T, the background-dependent path terms calculatedin step 806 of FIG. 8, and the candidate surface parameters (i.e.,directional emissivity, diffuse reflectance and Lambertian-equivalentsolar reflectance). In step 906, the candidate spectral radiancesignature and its uncertainty are compared with the target radiancesignature and its uncertainty calculated in step 808 of FIG. 8 tocalculate the delta signature and uncertainty.

In steps 908-912 of FIG. 9, a chi-square closeness-of-fit statistic,Jacobian matrix, and delta temperature δT are respectively calculatedfor the delta signature and its uncertainty. In step 914, adetermination is made as to whether the delta temperature δT is smallenough to indicate a match between the target radiance signature and thecandidate spectral radiance signature, thus indicating that thecandidate spectral radiance signature is, in fact, the best-fittingcandidate spectral radiance signature. If the delta temperature δT istoo large, the candidate spectral radiance signature is assumed not tobe the best-fitting candidate spectral radiance signature, thetemperature T is updated in step 916, and steps 904-914 are repeateduntil the best-fitting candidate spectral radiance signature isdetermined.

FIG. 10 illustrates a second method, namely, a hyperspectral method, forcalculating the best-fitting candidate spectral radiance signature ofstep 814. In FIG. 10, the Iterative Spectrally SmoothTemperature-Emissivity Separation (ISSTES) method described in C. C.Borel, “Iterative retrieval of surface emissivity and temperature for ahyperspectral sensor”, Proc. 1st JPL Workshop Remote Sensing of LandSurface Emissivity, published in 1997 by the Jet Propulsion Laboratory,is used to estimate the target temperature and the target spectralemissivity. In step 1002 of FIG. 10, the background-dependent path termsand target spectral radiance signature are provided as input to theISSTES method to calculate the target temperature. In step 1004, thecalculated target temperature and candidate surface parameters (i.e.,directional emissivity, diffuse reflectance and Lambertian-equivalentsolar reflectance) are then used to calculate a best-fitting candidatespectral radiance signature that corresponds to the target temperatureand candidate surface parameters.

Referring back to FIG. 8, in step 816, a computer processor may be usedto calculate the difference and error for the best-fitting candidatespectral radiance signature estimated in step 814. At step 818, theprocessor uses a threshold computed from the desired value for the P_(d)(probability of detection) to determine whether the best-fittingcandidate radiance of the target matches within a predetermined level ofconfidence of one of the one or more candidate spectral radiancesestimated at step 814. Candidate spectral radiances that match are usedto generate a list of possible target materials.

Polarimetric Reflective Processing

Referring to FIG. 11, imagery generated by radiation reflecting off thephysical target and sensed by a polarimetric sensor is referred to aspolarimetric reflective imagery. The polarimetric sensor typicallyincludes four orientations of polarization filters: 0°, 45°, 90°, and135°, although it could be done with only three polarization filters.The polarimetric sensor measures a Stokes vector in each spectral band.Polarimetric reflective imagery may be single band, but is preferablyeither multi-band or hyperspectral. Because light reflected off thetarget and sensed by the polarimetric sensor may be incoherent, thetarget signature may be measured in a Mueller matrix parameter spectrum.Error sources for polarimetric reflective imagery may includeatmospheric conditions, radiometric calibration, BRDF measurements,polarimetric calibration, and spectral calibration.

FIG. 11 illustrates an example of an implementation of the materialidentification process of FIG. 2 for polarimetric reflective imagery.For each image data set, the example in FIG. 11 utilizes the followinginputs: atmospheric information (i.e., information about the atmospherethrough which the radiation was sensed), candidate material surfaceproperties, uncertainty inputs, relative spectral response functions,and background and target radiance measurements.

Data regarding the atmosphere includes the same aerosol parameters asdescribed above in the non-polarimetric reflective process. For apolarimetric sensor, the relative spectral response functions typicallyconsist of four functions per band—one for each orientation of thepolarization filters (0°, 45°, 90°, and 135°). Additionally, thecandidate material surface properties consist of a Mueller matrix ineach band. The background and target radiance measurements consist ofStokes vectors in each band instead of the scalar measurements providedin the non-polarimetric reflective case. The uncertainty inputs includeradiance measurement uncertainty, candidate material Mueller matrixuncertainty, aerosol parameter uncertainty, and spectral calibrationuncertainty (for a hyperspectral case). The radiance measurementuncertainty comprises a covariance matrix for errors between the Stokesvector components in each band and, in the multispectral polarimetriccase, for the errors between the bands. The candidate material Muellermatrix uncertainty comprises a covariance matrix for errors between allentries in the Mueller matrix, instead of for the diffuse and solarreflectance in the non-polarimetric reflective case.

In step 1102, the aerosol parameters are provided as input to determinebackground-independent atmospheric correction terms. Thebackground-independent atmospheric correction terms are used inconjunction with the measured image radiance data to estimate backgroundreflectance as shown in step 1104, assuming that the background consistsof a non-polarizing material.

In step 1106, background-dependent atmospheric correction terms areestimated using the background reflectance data as well as the aerosolparameters. In step 1108, the process uses the input data to estimatethe BRDF atmosphere needed for estimating the Mueller matrix. In step1110, one or more candidate materials are selected from the BRDFdatabase 1111; and the process estimates the Mueller matrix, which isused with the background-dependent path terms to predict candidateStokes vectors and error bars in step 1112. The predicted radiance for acandidate material, in this example, consists of a Stokes vector in eachband (instead of the scalar values for the non-polarimetric reflectivecase described above). Calculation of these values involves the Stokesvectors describing the polarization state of the upwelling radiance, andoperations with Mueller matrices for the material and the Stokes vectorsdescribing the polarization state of the downwelling radiance.

In step 1114, the measured image radiance is used to estimate a targetStokes vector and its corresponding error bars. Finally, at step 1116, acomputer processor may be used to calculate the difference and error forthe predicted candidate Stokes vector and the estimated target Stokesvector. At step 1118, the processor uses a threshold computed from thedesired value for the probability of detection (P_(d)) to determinewhether the predicted candidate Stokes vector indicates a possiblematerial matching the target.

Polarimetric Emissive Imagery

Referring to FIG. 12, imagery generated by sensing the emissivity of thephysical target using a polarimetric sensor is referred to aspolarimetric emissive imagery. The polarimetric sensor typicallyincludes four orientations of polarization filters: 0°, 45°, 90°, and135°, and measures a Stokes vector in each spectral band. Polarimetricemissive imagery may be single band, but is preferably either multi-bandor hyperspectral. Error sources for polarimetric reflective imagery mayinclude atmospheric conditions, radiometric calibration, BRDFmeasurements, polarimetric calibration, and spectral calibration.

FIG. 12 illustrates an example of processing polarimetric emissiveimagery in accordance with the material identification process of FIG.2. For each image data set, the example in FIG. 12 utilizes thefollowing inputs: relative spectral response functions similar to thoseprovided in the polarimetric reflective case, atmospheric informationsimilar to that provided in the non-polarimetric emissive case,background and target radiance measurements similar to those provided inthe polarimetric reflective case, candidate material surface properties,and uncertainty inputs.

The candidate material surface properties include a Mueller matrix ineach spectral band such as that in the polarimetric reflective case, andpolarized directional emissivity. Uncertainty inputs include radiancemeasurement uncertainty such as that provided in the polarimetricreflective case, candidate material property uncertainty comprising acovariance matrix for the errors between all entries in the Muellermatrix, a covariance matrix for the emissivity values in each band, anda covariance matrix for errors between the Mueller matrix entries andthe emissivity values, atmospheric parameter uncertainty such as thatprovided in the non-polarimetric case, and spectral calibrationuncertainty for the hyperspectral case.

In step 1202 of FIG. 12, a background-independent atmosphericcorrection, comprised of background-independent path terms, is estimatedusing the input data. The background-independent atmospheric path termsare combined with the image radiance data to estimate backgroundreflectance and emissivity of the target as shown at step 1204, assumingthat the background consists of a non-polarizing material.

At step 1206, the input data, background reflectance and emissivity, andbackground-independent path terms are combined to estimate abackground-dependent atmospheric correction comprised ofbackground-dependent path terms. At step 1208, the image radiance datais used to estimate the target Stokes vector and the error correspondingto the target Stokes vector.

At step 1210, the process uses the input data to estimate the BRDFatmosphere needed for estimating the Mueller matrix and emissivity. Atstep 1212, one or more candidate materials are selected from BRDFdatabase 1213; then, for each candidate material, a Mueller matrix andemissivity are estimated. The predicted radiance for a candidatematerial comprises a Stokes vector in each band, wherein calculation ofthese values involves operations similar to those used for thepolarimetric reflective case, plus the calculation of the thermalemission contributions from the atmosphere and the target. In step 1214,the best-fitting candidate Stokes vector and its error are determined bycomparing the candidate Stokes vector to the estimated target Stokesvector.

At step 1216, a computer processor may be used to calculate thedifference and error for the best-fitting candidate Stokes vector andthe estimated target Stokes vector. At step 1218, the processor uses athreshold computed from the desired value for the probability ofdetection (P_(d)) to determine whether the best-fitting candidate Stokesvector indicates a possible material matching the target.

Synthetic Aperture Radar (SAR) Imagery Processing

Referring to FIG. 13, SAR imagery is generated by sensing radiance ofthe physical target in the microwave region of the spectrum. For SARimagery, two configurations are possible: mono-static SAR andmulti-static SAR. For a mono-static SAR configuration, at least oneimage is taken near the normal orientation of the target surface, and atleast one image is taken far from the normal orientation. For amulti-static SAR configuration, at least one image is taken near thespecular direction, and at least one image is taken far from thespecular direction. The target signature for SAR imagery may be measuredin terms of the radar cross-section spectrum. Accordingly, error sourcesfor SAR imagery may include atmospheric conditions, radiometriccalibration, and BRDF measurements.

FIG. 13 illustrates processing of SAR imagery in the radiance domain.The inputs specific to this example include the following: targetradiance measurements consisting of Stokes vectors in each band,candidate material surface properties consisting of the radarcross-section for a material with respect to the polarization state ofthe incident and reflected waves, and uncertainty inputs. Theuncertainty inputs include radiance measurement uncertainty andcandidate material property uncertainty consisting of a covariancematrix for the errors between all entries in the Mueller matrix such asfor the polarimetric reflective case. A primary difference between SARand the other modalities is the presence of specular noise resultingfrom coherent interference during image processing.

In step 1302 of FIG. 13, a radar cross-section is estimated using one ormore candidate materials selected from the BRDF database 1303. Thepredicted radiance for a candidate material, in this example, consistsof a Stokes vector in each band. These values are calculated usingmicrowave diffraction theory. In step 1304, a target Stokes vector andits error are estimated using the measured image radiance data. In step1306, the candidate Stokes vector and its error are predicted using theradar cross-section data estimated in step 1302. As shown in step 1308,by calculating the difference in the signatures of the candidate Stokesvector and the estimated target Stokes vector, it may be determined thatthe material selected from the BRDF database may be added to the list ofpossible target materials, as shown at step 1310.

The foregoing description is of exemplary and preferred embodimentsemploying at least in part certain teachings of the invention. Theinvention, as defined by the appended claims, is not limited to thedescribed embodiments. Alterations and modifications to the disclosedembodiments may be made without departing from the invention as setforth in the amended claims. The meaning of the terms used in thisspecification are, unless expressly stated otherwise, intended to haveordinary and customary meaning and are not intended to be limited to thedetails of the illustrated structures or the disclosed embodiments.

What is claimed is:
 1. A method performed by processing circuitry andmemory configured for determining if a physical target is formed from acandidate material, comprising: receiving a first image of the physicaltarget acquire first acquisition angle; delineating the physical targetwithin the first image as a first area; extracting a first spectralradiance value from the first area; receiving a second image of thephysical target acquired at a second acquisition angle different fromthe first acquisition angle; delineating the physical target within thesecond image as a second area; extracting a second spectral radiancevalue from the second area; evaluating an estimated target spectralsignature for the physical target based at least in part on the firstand second spectral radiance values; calculating a predicted candidatespectral signature for the candidate material based at least in part onthe first and second acquisition angles; calculating a correlationstatistic that correlates the predicted candidate spectral signature tothe estimate target spectral signature; comparing the correlationstatistic to a predetermined probability threshold; and determining fromthe comparison that the physical target is formed from the candidatematerial.
 2. The method of claim 1, wherein the first acquisition angleis within a forward scattering lobe of a radiation being reflected fromthe target from a known illumination source, and the second acquisitionangle is outside the forward scattering lobe.
 3. The method. of claim 1,wherein at least one of the first or second images is acquired using anon-polarimetric reflectance imaging modality.
 4. The method of claim 1,wherein calculating the spectral signature comprises obtaining surfacereflection parameters for the candidate material using a bidirectionalreflectance distribution function (BRDF).
 5. The method of claim 1,further comprising determining errors in the predicted candidatespectral signature based one or more uncertainties.
 6. The method ofclaim 5, wherein calculating the predicted candidate spectral signaturefurther comprises use of the one or more surface reflection parametersobtained from a predetermined bidirectimial reflectance distributionfunction (BRDF) for the candidate material, and wherein the one or moreuncertainties include one or more uncertainties in the one or moresurface reflection parameters for the candidate material.
 7. Anon-transitory computer readable media storing instructions that, whenexecuted by a computing system, cause the computing system to perform amethod for determining if a physical target is formed from a candidatematerial, the method comprising: receiving a first image of the physicaltarget acquired at a first acquisition angle; delineating the physicaltarget within the first image as a first area; extracting a firstspectral radiance value from the first area; receiving a second image ofthe physical target acquired at a second acquisition angle differentfrom the first acquisition angle; delineating the physical target withinthe second image as a second area; extracting a second spectral radiancevalue from the second area; evaluating an estimated target spectralsignature for the physical target based at least in part on the firstand second spectral radiance values; calculating a predicted candidatespectral signature for the candidate material based at least in part onthe first and second acquisition angles; calculating a correlationstatistic that correlates the predicted candidate spectral signature tothe estimated target spectral signature; comparing the correlationstatistic to a predetermined probability threshold; and determining fromthe comparison that the physical target is funned from the candidatematerial.
 8. The non-transitory computer readable media of claim 7,wherein the first acquisition angle is within a forward scattering lobeof a radiation being reflected from the, target from a knownillumination source, and the second acquisition angle is outside theforward scattering lobe.
 9. The non-transitou computer readable media ofclaim 7, wherein at least one of the first or second images is acquiredusing a non-polarimetric reflectance imaging modality.
 10. Thenon-transitory computer readable media of claim 7, wherein calculatingthe spectral signature comprises obtaining surface reflection parametersfor the candidate material using a bidirectional reflectancedistribution function (BRDF).
 11. The non-transitory computer readablemedia of claim 7, wherein the method further comprises determiningerrors in the predicted candidate spectral signature based one or moreuncertainties.
 12. The non-transiton computer readable media of claim11, wherein calculating the predicted candidate spectral signaturefurther comprises use of the one or more surface reflection parametersobtained from a predetermined bidirectional reflectance distributionfunction (BRDF) for the candidate material, and wherein the one or moreuncertainties include one or more uncertainties in the one or moresurface reflection parameters for the candidate material.
 13. Anapparatus for determining if a physical target is formed from acandidate material, comprising: a specially programmed computing system,the specially programmed computing system including a processor, adisplay, and memory configured to store instructions that are executableon the processor, the instructions configured to: receive a first imageof the physical target acquired at a first acquisition angle; delineatethe physical target within the first image as a first area; extract afirst spectral radiance value from the first area; receive a secondimage of the physical target acquired at a second acquisition angledifferent from the first the first acquisition angle; delineate thephysical target within the second image as a second area; extract asecond spectral radiance value from the second area; evaluate anestimated target spectral signature for the physical target based atleast in part on the first and second spectral radiance values;communicate with a bidirectional reflectance distribution function(BRDF) library; calculate a predicted candidate spectral signature forthe candidate material based at least in part on the first and secondacquisition angles and based at least in part on the BRDF library;calculate a correlation statistic that correlates the predictedcandidate spectral signature to the estimated target spectral signature;compare the correlation statistic to a predetermined probabilitythreshold; and determining from the comparison that the physical targetis formed from the candidate material.
 14. The apparatus of claim 13,wherein the first acquisition angle is within a forward scattering lobeof a radiation being reflected from the target from a known illuminationsource, and the second acquisition angle is outside the forwardscattering lobe.
 15. The apparatus of claim 13, wherein at least one ofthe first or second images is acquired using a non-polarimetricreflectance imaging modality.
 16. The apparatus of claims 13, whereinthe instructions are further configured to determine errors in thepredicted candidate spectral signature based one or more uncertainties.17. The apparatus of claim 16, wherein calculating the predictedcandidate spectral signature uses of the one or more surface reflectionparameters obtained from a predetermined bidirectional reflectancedistribution function (BRDF fur the candidate material, and wherein theone or more uncertainties include one or more uncertainties in the oneor more surface reflection parameters for the candidate material.